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SUMMARY 


A new numerical framework for solving conservation laws is being developed. This new 
approach differs substantially from the well established methods, i.e., finite difference, finite 
volume, finite element, and spectral methods, in both concept and methodology. It employs: 

a. a nontraditional formulation of the conservation laws in which space and time are unified 
and treated on the same footing; and 

b. a nontraditional use of discrete variables such that numerical marching can be carried out 
by using a set of relations that represents both local ami global flux conservation. 

To be specific, we consider a conservation law that governs the convection and diffusion of a 
physical variable in a 1-D space. Let (i) x be the spatial coordinate, (ii) cq > 0 be a conversion 
constant with the dimension of velocity, and (iii) t be the product of c 0 and the temporal 
coordinate. By definition, x and t have the same dimension. As a result, x\=x and x 2 — f rnay be 
considered as the coordinates of a two-dimensional Euclidean space E 2 (also referred to as a 
space-time). Let u(x,t) be a scalar function of x and t. Let a be a dimensionless constant and p 
(> 0) be a constant with the dimension of length. Then the conservation law may be expressed 
as 



where (i) S (V) is the boundary of an arbitrary volume V in E 2 , (ii) 


( 0 . 1 ) 


def 


du 


5 <— *£■“> 


( 0 . 2 ) 


is a current density vector in E 2 , and (iii) ds = dolt with do and it, respectively, being the area 
and the outward unit normal of a surface element on S (V). By applying the divergence theorem 
in£ 2 . Eq. (0-1) implies the unsteady convection-diffiision equation, i.e.. 


du du d 2 u _ 

-T- + a — - p — r = 0 

dt dx dx 2 


(0.3) 


Let E 2 be divided into nonoverlapping regions referred to as conservation elements (see Figs. 
2.1(a) and 2.1(b)). A conservation element and its interior, respectively, may be denoted by 
CE(j,n) and CE"(j,n) where j and n, respectively, are the spatial and temporal indices. For 
(x,t) e CE"(J,n), u(x,t ) will be approximated by 

u(x,t) a j(x-x") + P"(r-r") + yj (0.4) 

where a", P" and are constants in CE "(j,n), and (x",t n ) are the coordinates of the center of 
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(0.5) 


CE(j\n). For (x f t) e CE"(J f n ), h (x,t) will be approximated by 

=? {an(x,t)-]X. ^ , n (x, t ) ) 

OX 

Furthermore, the conservation law Eq. (0. 1) will be approximated by 

W^ = 0 (a6) 

where V is the union of any combination of conservation elements. Since 71* is not defined on 
the above surface integration, by definition, is to be carried out over a surface which is in 
the interior of V and immediately adjacent to S (V). 

Because and li(x,t) are continuous in the interior of a conservation element but may be 
discontinuous across an interface separating two neighboring conservation elements, a 
conservation element is also a solution element in the current scheme. As will be shown, 
generally a conservation element is not necessarily a solution element and vice versa. 

Let V= CE (j,n). Then Eqs. (0.4) - (0.6) and the divergence theorem imply that P" = -a a". 
As a result, Eq. (0.4) can be simplified as 

u(x,t) = a] [(x-Xj)-a(t-t n )] + i] , (x,t) e CE"(j,n) (0.7) 

Thus, for ( x,t ) e CE u(x,t ) is determined by the parameters y y and a". As will be shown, 
by repeatedly applying Eq. (0.6) with V being the union of two neighboring conservation 
elements, any pair of y ” and a] can be determined in tenms of and a”, j = 0, ±1, ±2, • • • . The 
values of and a], in turn, can be determined by the initial condition. 

Because 


u(x],t n ) = y 1 (0.8) 

and 

du 

-jjj- = a j (x,t)e CE"(j,n) (0.9) 

Yj and a ", respectively, may be considered as the numerical counterparts of u(x],t n ) and 
u x (xj,t n ). In other words, both u and its spatial derivative at (Xj,t n ) are computed by the current 
scheme. 

In the current paper, we also 

a. explore the concept of a dynamic space-time mesh (the conservation elements are 
embedded in this mesh) and the need for a unified treatment of physical variables and mesh 
parameters; 
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b. study the stability, dissipation and dispersion of the current scheme by using a rigorous 
Fourier analysis; 

c. develop a new error analysis technique that enables us to predict and interpret the 
numerical errors of the current and other classical schemes; 

d. study the consistency and truncation error of the current scheme; and 

e. compare the errors of the numerical solutions generated by the current scheme and other 
classical schemes. 

The key results obtained from the above study are: 

a. It is demonstrated that (i) stability and accuracy can be improved, and (ii) dissipation and 
dispersion can be reduced, if the space-time mesh is allowed to evolve with the physical 
variable such that the local convective motion of physical variables relative to the moving 
mesh is kept to a minimum. Because the appearance of wiggles near a discontinuity is a 
result of numerical dispersion, these wiggles can also be reduced by reducing this relative 
convective motion. 

b. It is shown that there is a remarkable similarity between the forms of the amplification 
factors of the Leapfrog/ DuF ort-Frankel and the current schemes. As a result of this 
similarity, the stability condition of the current scheme, as in the case of the 
Leapfrog/DuFort-Frankel scheme, is essentially the CFL condition and thus independent of 
the viscosity p. Therefore, the current scheme is unconditionally stable in the case of pure 
diffusion. Also, as in the case of the Leapfrog/DuFort-Frankel scheme, the current scheme 
has no numerical diffusion in the absence of viscosity. Note that the stability condition of a 
classical explicit scheme for solving Eq. (0.3), e.g., the MacCormack scheme, generally is 
more restrictive than the CFL condition. In the case where the mesh Reynolds number <e 
1, the stability bound for the time-step size At is more or less proportional to (Ax) 2 . In 
contrast, the same bound will still be determined by the CFL condition and therefore is 
proportional to Ax if the current scheme is used. The advantage of the current scheme in 
the allowable time-step size grows as Ax 0. This advantage becomes particularly 
important when the current scheme is used in a steady-state calculation. 

c. It is shown theoretically that the current scheme is more accurate than the 
Leapfrog/DuFort-Frankel scheme by one order (in a sense to be defined in the paper) in 
both initial-value specification and the marching scheme. It is also shown theoretically that 
the current scheme is substantially more accurate than the MacCormack scheme in spite of 
their almost identical operation counts. Its advantage ranges from a factor of four for the 
case of pure convection to several orders of magnitude if diffusion is dominant and a 
theoretically-determined optimal At is used. 
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d. It is shown that the consistency of the current scheme, as in the case of the 
Leapfrog/DuFort-Frankel scheme, requires that At /Ax -> 0 as At, Ax 0. This fact 
contrasts sharply with most other explicit schemes, e.g., the MacCormack scheme, that 
have no such requirement for consistency. However, by using Lax’s equivalence theorem 
and a necessary condition for convergence, it is shown that, for these explicit schemes, this 
requirement must manifest itself as a part of stability conditions. As a matter of fact, it is 
shown that the truncation errors of the Leapfrog/DuFort-Frankel. the MacCormack, and 
the current schemes are all second order in Ax if stability is taken into consideration. 

e. It is shown numerically that the current scheme is far superior to the Leapfrog/DuFort- 
Frankel scheme in accuracy, and has a substantial advantage over the MacCormack scheme 
in both accuracy and stability. 



1. INTRODUCTION 


This paper is the first of a series of papers that describe a new framework for solving 
conservation laws. This new approach differs substantially from the well established methods, 
i.e., finite difference, finite volume, finite element, and spectral methods, in both concept and 
methodology. The focus of the current numerical simulation is entirely on the integral forms of 
the conservation laws. Little or no attempt is made to simulate the differential forms which are 
valid only when the dynamical variables are well behaved. As a result, this new framework has 
the potential to provide more accurate simulation of the physical phenomena in which the 
dynamical variables may not vary smoothly. 

Specifically, the explicit scheme to be presented in this paper employs: 

a. a nontraditional formulation of the conservation laws in which space and time are unified 
and treated on the same footing; and 

b. a nontraditional use of discrete variables such that numerical marching can be carried out 
by using a set of relations that represents both local and global flux conservation. 

As a preliminary, this paper will begin with a discussion on the conservation laws. For 
simplicity, we consider the conservation law that governs the convection of a physical variable in 
a 1-D space. Let (i) x be the spatial coordinate, (ii) c 0 > 0 be a conversion constant with the 
dimension of velocity, and (iii) t be the product of cq and the temporal coordinate. By definition, 
x and t have the same dimension. As a result, x x =x and x 2 = t may be considered as the 
coordinates of a two-dimensional Euclidean space £ 2 [p.161, 1]. Note that the scalar product of 
any two vectors in E 2 is defined in [1] as a part of the definition of £ 2 - Let u(x,t) be a scalar 
function of x and t. Let a be a dimensionless constant. Then the conservation law can be 
expressed as 



where (i) S (V) is the boundary of an arbitrary volume V in £ 2 , (H) 

li S' ( au , u) (1.2) 

is a current density vector in £ 2 , and (iii) ds = dolt with do and it, respectively, being the area 
and the outward unit normal of a surface element on 5 (V). Note that an n-dimensional Euclidean 
space E n (n > 2) may be referred to as a space-time if one of its coordinates is temporal in nature 
while others are spatial in nature. Also h • ds may be referred to as the flux of h leaving the 
volume V through the surface element ds. With the aid of the divergence theorem and the fact 
that rf li= d(au) /dx + du /5r, one may obtain the differential form of Eq. (1. 1), i.e., 
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du du 

— + a ^r- 
dt dx 


(1.3) 


= 0 


The current unified and equal treatment of space and time is in sharp contrast to the 
traditional approach in which space and time are divided and treated separately. The following 
comments are made to clarify the differences between these two approaches: 

a. Geometric figures referred to in the traditional approach, such as rectangles and triangles, 
usually are objects in space (Note: By definition, "space" is the coordinate hyperplane with 
time = 0. See p.263 in [2]). Contrarily, geometric figures referred to in the current paper, 
unless specified otherwise, are objects in space-time. Note that, in this paper, a geometric 
figure, such as a rectangle, implies both its boundary and interior. 

b. In a space-time E 2 , the volume V in Eq. (1.1) is traditionally taken to be a rectangle with 
its edges being aligned with either the x- axis or the t- axis. With this choice, the integral 
on the left side of Eq. (1.1) can be divided into four parts, each of which involves only 
integration in time or space. Contrarily, in the current approach, the volume V can be a 
geometric figure of any shape and thus the surface integration over the boundary of V may 
involve both space and time simultaneously. 


c. 


In a space-time E n with n> 3, we also can consider a conservation law with the form of Eq. 
(1.1). For example, the mass conservation law in a space-time E 3 can be expressed in the 
form of Eq. (1.1) with the coordinates x i , x 2 and x 3 , and the vector/i* defined by 





(1.4) 


and 


h = (— p , ~~p , p) 
Co C 0 


(1.5) 


Here (i) y is a spatial coordinate, (ii) p is the mass density, and (Hi) v x and v y are the 
velocity components in the x - and y- directions, respectively. In the traditional approach, 
the volume V in Eq. (1.1) is taken to be a cylinder in £3 like that depicted in Fig. 1 . 1 . 
Assuming that (i) each of the two ends of this cylinder has a constant value of t, and (ii) the 
generators of its side surface point in the t- direction, then Eq. (1.1) implies that 



(i+M) 



1+&1 


+ — J dt'f p^- ds ± = 0 
C 0 1 , 


( 1 . 6 ) 


where (i) V x is the projection of the cylinder on the x-y plane, (ii) dv x is a volume 
element in V x , (iii) Sx(V / _ L ) is the boundary of V x , (iv) ds x is a surface element on 
S.l(F-l), ^ ( v )^ d - ( y*, v y , 0 ) is a vector lying on the x-y plane. Traditionally, V x is 
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referred to as the "control volume" and Eq. (1.6) is known as the integral form of the two- 
dimensional time-dependent mass conservation law. Note that, by definition, t I cq is the 
ordinary temporal coordinate. Since the control volume V± is an object in space, the first 
two integrals in the above equation involve only integration in space while the last integral 
represents a combined operation in which a surface integration in space is followed by an 
integration in time. This is another example in which space and time are divided and 
treated separately. With the division, the above equation has a simple interpretation, i.e., 
the increase of the mass in the control volume V± during a time interval At /cq is equal to 
the total mass entering V j_ through its boundary Si(V_l) during the same time interval. 
However, the division of space and time is achieved at the expense of limiting the choice of 
the volume V to a cylinder in space-time like that depicted in Fig. 1.1. Note that pl^ a 
vector in space, is commonly known as the mass current density vector. Also pT*- ds± is 
referred to as the flux of p”v^ leaving the control volume through the surface element 
ds±. These definitions involve only vectors in space. On the contrary, the current density 
vector h is a vector in space-time while the flux li ds is the inner product of two vectors in 
space-time. 


The above remarks make it clear that a greater flexibility in choosing the volume V is allowed 
in the current formulation of conservation laws than that allowed in the traditional approach. As 
will be shown, the use of this flexibility is an integral part of the current numerical framework. 

Next the classical Lax-Wendroff scheme will be discussed using the uniform mesh depicted 
in Fig. 1.2(a). This discussion is presented such that readers may understand the stream of 
thoughts that leads to the development of the current framework. 


The Lax-Wendroff scheme for solving Eqs. (1.1) and (1.2) consists of two distinctly different 
marching steps. In the first step, the variables Uj$ at the time level (n+Vi) are evaluated in 
terms of the variables u" at the time level n, i.e., 


U j+'A 



(1 +v)u" + (1-V 



(1.7) 


where 

v a At / Ax (1.8) 

is the Courant number. The derivation of Eq. (1.7) may be explained using Fig. 1.2(a). Point Q 
is at the time level « and on the same characteristic line with the mesh point P. The value of u at 
P, i.e., u"t $ is evaluated in terms of those at the mesh points R and S, i.e., uj+i and u], by 
assuming (i) u is constant along a characteristic line, and (ii) the value of u at Q is the linear 
interpolation of those at R and S. 
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In the second step, the variables u] +l at time level n + 1 are determined in terms of the 
variables at time levels n and n+Vi by using a conservation relation. Specifically, one assumes 
that Uj , Wy+v 4 1 Uj and , respectively, represent the average values of u on the line segments 

BC, CD, DA and AB. Then an application of Eq. (1.1) with V being the rectangle ABCD implies 
that 


u ] +1 Ax - u" Ax + a u n j$A At - a u At = 0 (1.9) 

Eq. (1.9) is equivalent to 

«; +1 - U ]%) ( 1 . 10 ) 

This is the relation used in the second marching step. Substituting Eq. (1.7) into Eq. (1.10), one 
obtains a difference form in which the variables at time level n+l are expressed directly in terms 
of the variables at time level n, i.e., 


«r' - ♦ (i 


ci.ii) 


To serve as the starting point of the current development, the conservation relation (1.9) will 
be cast into a form similar to Eq. (1 . 1). To proceed, let (see Fig. 1.2(b)) 

f u] if (x,t) e 0" IBGC 

u(x,t) \ (1.12) 

[ u n jt£ if (x,t) e 0" DICF 

where 0" IBGC denotes the interior of the rhombus 0 IBGC, and so on. Similarly, u(x,t) can be 
defined for any (x,t) which is in the interior of other similar rhombuses. Since u(x,t) is 
continuous in the interior of each of these rhombuses but may be discontinuous across an 
interface separating two neighboring rhombuses, such a rhombus will be referred to as a solution 
element. In terms of u(x,t), the vector function ]i(x,t) is defined by 

?(*.0 - (au(x,t) , «(x,0) (1.13) 

With the aid of Eqs. (1.12) and (1.13), Eq. (1.9) can be expressed as 

F ) 

9 h - ds = 0 (l 141 

J S(OABCD) - 

i&., the total flux leaving the rectangle □ ABCD vanishes if /ns the flux density vector. Note that 
h is not defined at the vertices A, B, C and D. However, contributions to the above integral from 
these isolated points are zero no matter what /fare assigned to them. As a result, they may simply 
be excluded from the above surface integration. Since Eq. (1.9) applies for any 
j = 0* and « =0,1,2, ■■■, Eq. (1.14) is valid if □ ABCD is replaced by any similar 

rectangle like □ JADK or □ DCML shown in Fig. 1 .2(b). For this reason, each of these rectangles 
will be referred to as a conservation element for the Lax-Wendroff scheme. Note that, excluding 
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its two end points, an interface separating two neighboring conservation elements is located in the 
interior of a solution element. As a result, the vector function It is continuous across such an 
interface. This coupled with Eq. (1.14) implies that the total flux leaving any volume V which is 
the union of any combination of conservation elements must also vanish, i.e.. 



For example, V can be the L-shaped figure formed by □ ABCD, □ JADK and DDCML. Eq. 
(1.15) is in a form similar to Eq. (1.1). However, it is equivalent to Eq. (1.9) and thus represents 
only one of two marching steps that form the Lax-Wendrojf scheme. 

At this juncture, we emphasize that both solution elements and conservation elements are 
domains in space-time. Contrarily, elements in the finite element method are domains in space 
only. 

In its earliest form, the current scheme may be considered as a modification of the Lax- 
Wendroff scheme in which all the marching steps are derived from a single conservation 
relation. The modifications begin with the assumption that u(x,t ) is approximated by 

u(x,t) = a ](x-x]) + P ](t-t n ) + fj if (x,t) g 0"IBGC (1.16) 

where (i) (x],t n ) are the coordinates of the center of OlBGC depicted in Fig. 1.2(b), and (ii) a], 
P" and i] are considered constants in 0" IBGC. Note that x* is only a function of j if a 
"stationary" space-time mesh, e.g., a mesh shown in Fig. 1.2(b), is considered. However, in 
Section 2, it becomes a function of both j and n when a "moving" mesh is introduced. Also note 
that 

u(x],t n ) = y ] (1.17a) 

and 

|j = a j and = P? if (x,t)e 0" IBGC (1.17b) 

i.e., yj is the value of u at the center of 0 IBGC while a" and P", respectively, are the spatial and 
temporal derivatives of n in 0" IBGC. 

Hereafter, unless specified otherwise, an equation like Eq. (1.16) is assumed to be valid for 
any ( j,n ) with either (i) n = 0, 1, 2, •■•,) = 0, ±1, ±2, • • • , or (ii) n = 1/2, 3/2, 5/2, 
j =±1/2, ±3/2, ±5/2, • • • . Thus the rhombuses referred to earlier are also the solution elements 
in the current method. 

In the current method, it is also assumed that 
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(U8) 


(6 Ji- ds = 0 

where 

(au(x,t ) , u(x,t)) (1.19) 

is the flux density vector, and V is the union of any combination of the rhombuses referred to 
earlier. Since h is not defined on S (V), the above surface integration, by definition, is to be 
carried out over a surface which is in the interior of V and immediately adjacent to S(V). A 
necessary condition for the conservation relation Eq. (1.18) is 



where 0 is any one of the rhombuses referred to earlier. Thus these rhombuses are also the 
conservation elements in the current scheme. This is different from the Lax-Wendroff scheme in 
which a conservation element is a rectangle like □ ABCD depicted in Fig. 1.2(b). 

Another necessary condition for Eq. (1.18) is the requirement that the net flux of /^entering an 
interface separating two neighboring conservation elements (i.e., the rhombuses) must vanish. 
This may be seen by applying Eq. (1.18) separately to two neighboring rhombuses and then to the 
union of them. Obviously the local flux conservation relations at each interface and within each 
conservation element (i.e., Eq. (1.20)) are equivalent to the conservation relation Eq. (1.18). In 
the next section, the current marching procedure will be constructed by using the local 
conservation relations. 

This completes the description of the basic concepts behind the current development. In this 
first paper, these concepts will be used to construct a numerical scheme for solving an unsteady 
1-D constant-coefficient convection-diffusion model equation over a uniform constant-velocity 
moving mesh. The model equation and the mesh used are simple enough such that the important 
properties of the resulting scheme may be studied analytically. Yet they are complicated enough 
that the information gained and the techniques developed in the current study may provide a solid 
base for the development of new schemes for solving nonlinear conservation laws in higher 
dimension. Note that it has been shown empirically that the local behaviors of a nonlinear 
variable-mesh scheme may be studied by using a local analysis (such as the von Neumann 
analysis) in which the dynamic coefficients and geometric parameters are frozen at their local 
values. In the same spirit, the current analysis is intended to serve as a guide for the local 
analysis of the more complicated schemes to be developed later. 

The remainder of this paper is briefly described as follows: In Section 2, we construct the 
current scheme without using several questionable assumptions commonly made in the 
construction of an explicit, time-accurate, conservative scheme. We also point out several 
fundamental differences that separate the current scheme from the traditional schemes. One of 
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them is the fact that the current scheme is locally implicit while globally explicit. It is also 
explained how these differences will result in greater stability and accuracy for the current 
scheme. 

In Section 3, we explore the concept of a dynamic space-time mesh and the need for a unified 
treatment of the physical variables and mesh parameters. Specifically, it is demonstrated that the 
stability and accuracy of a numerical calculation may be improved if the space-time mesh is 
allowed to evolve with the physical variables such that the local convective motion of the 
physical variables relative to the moving mesh is kept to a minimum. In the meantime, a 
parameter defined in Section 2 is interpreted as the Courant number for a moving mesh. 

In Section 4, the stability, dissipation and dispersion of the current scheme are studied using a 
rigorous Fourier analysis. It is shown that there is a remarkable similarity between the forms of 
the amplification factors of the Leapfrog/DuFort-Frankel [p. 161, 3] and the current schemes. 
Note that, hereafter, the former will be referred to as the L/D-F scheme. As a result of this 
similarity, the stability condition of the current scheme, as in the case of the L/D-F scheme, is 
essentially the CFL condition and thus independent of the viscosity coefficient p. Therefore, the 
current scheme is unconditionally stable in the case of pure diffusion. Also, as in the case of the 
L/D-F scheme, the current scheme has no numerical diffusion in the absence of viscosity. Note 
that the stability condition of a classical explicit scheme for solving Eq. (2.2), e.g., the 
MacCormack scheme [p.163, 3], generally is more restrictive than the CFL condition (see Fig. 
4.1). In the case that the mesh Reynolds number «*: 1, the stability bound for At is more or less 
proportional to (Ax) 2 . In contrast, the same bound will still be determined by the CFL condition 
and therefore is proportional to Ax if the current scheme is used. The advantage of the current 
scheme in the allowable time-step size grows as Ax — » 0. This may become particularly 
important when the current scheme is used in a steady-state calculation. 

In Section 5, assuming smooth and periodic initial data, an error analysis technique is 
developed using the discrete Fourier analysis formulated in Section 4. The main achievement in 
this development is the derivation of a simple formula for predicting the numerical errors of the 
current scheme. This formula contains a principal part and a spurious part. The principal part 
grows linearly with the time-step number n while the spurious part is independent of n. Thus the 
principal part will become dominant as n increases. Furthermore, it will be shown that this 
error-prediction formula is valid up to any n as long as the numerical solution is still accurate up 
to this n. Similar error-prediction formulae are also given for the L/D-F and the MacCormack 
schemes. The prediction formula for the L/D-F scheme also contains a principal part and a 
spurious part while that for the MacCormack scheme contains only the principal part. By using 
these formulae, it will be shown that the current scheme is more accurate than the L/D-F scheme 
by one order (in a sense to be defined later) in both initial-value specification and the marching 
scheme. These formulae may also be used to show that the current scheme is substantially more 
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accurate than the MacCormack scheme. This section is concluded by showing that the operation 
counts for the current scheme and the MacCormack scheme are almost identical. 

In Section 6, it is shown that the consistency of the current scheme, as in the case of the L/D - 
F scheme, requires that At/Ax ► 0 as At, Ax — > 0. This contrasts sharply with most other explicit 
schemes, e.g., the MacCormack scheme, which have no such requirement for consistency. 
However, by using Lax’s equivalence theorem [p.45, 4] and a necessary condition for 
convergence, it is shown that, for these explicit schemes, this requirement must manifest itself as 
a part of the stability conditions. As a matter of fact, it is shown that the truncation errors of the 
MacCormack, the L/D-F, and the current schemes are all second order in Ax if stability is taken 
into consideration. 

In Section 7, numerical solutions generated by the MacCormack, the L/D-F, and the current 
schemes are compared with the corresponding analytical solutions for different values of physical 
coefficients, mesh parameters and total running time. These comparisons show that the current 
scheme is far superior than the L/D-F scheme in accuracy, and has a substantial advantage over 
the MacCormack scheme in both accuracy and stability. Moreover, they confirm many of the 
theoretical predictions made earlier in this paper. 

Finally, odds and ends are dealt with in Section 8. They include discussions on boundary- 
value specification, conservation elements of other geometric shapes, and a possible extension of 
the current scheme to a space-time of higher dimension. 


2. MARCHING PROCEDURES 


In Section 1, for simiplicity, we consider only pure convection. Thus the flux density vector 
h—{au,u). Hereafter both convection and diffusion will be considered. As a result, Eq. (1.2) is 
replaced by 

T* def , du 

h ^ (au - ii — , u ) (2.1) 


where p > 0 is a constant with the dimension of length. We will continue to assume Eq. (1.1) and 
the related assumptions. Note that the unsteady convection-diffusion equation 


du du 

— + a — 
dt dx 



( 2 . 2 ) 


follows from Eq. (1.1) and (2.1) if u is well-behaved. 


For reasons to be explained later, we will also consider a moving mesh shown in Fig. 2.1(a). 
By "moving mesh", we mean a space-time mesh such that the coordinate x may vary along a j 
mesh line, i.e., a mesh line with a constant value of the index j. In Fig. 2.1(a), b is a constant and 
dx/dt = b along any j mesh line. In other words, a particle with a space -time trajectory 
coinciding with a j mesh line has a constant velocity b. For this reason, b may be referred to as 
the velocity of the moving mesh. The moving mesh is reduced to a stationary mesh if b = 0. Let 
the origin of the coordinate system coincide with the mesh point with j = n= 0. Then the 
coordinates x and t for a mesh point (J,n) are given by 

x = x] =? jAx + n bAt and t = t n d - nAt (2.3) 

In the current method, a conservation element is a parallelogram like that depicted in Fig. 
2.1(b). It is also a solution element. Hereafter, no distinction will be made between a 
conservation element and a solution element. A conservation element which is centered at 
(*",/") will be denoted by CE(j,n). Its interior will be denoted by CE"(j,n). 

To construct the marching procedure, u(x,t ) is assumed to be in the form defined by Eq. 
(1.16) for (x,t) e CE" (/,«). Moreover, we assume the conservation relation Eq. (1.18) with 

li{x, t) ( uu (x,t) — p , u (x, t ) ) (2.4) 

OX 

where V is the union of any combination of conservation elements. Obviously, this assumption is 
again equivalent to Eq. (1.20) and the interface flux conservation condition referred to in Section 
1. The only modification required is the generalization of conservation elements from rhombuses 
to parallelograms. 
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With the above assumptions, one has 


1? ^ , dll X , du „ an 

? h j -) +lr .,o; + pj 


(2.5) 


Thus the diffusion term in h does not contribute to V • ^ However, the diffusion term does play a 
role in the flux balance across an interface. With the aid of Eq. (2.5) and the divergence theorem, 
the generalized form of Eq. (1.20) implies that aa] + P" =0. Asa result, Eq. (1.16) implies that 

U(x,t) = a.j[(x-Xj)-a(t-t n )] + y] if (x,r)e CE"(J,n) (2.6) 

As a preliminary to the application of the interface flux conservation relation, next we 
consider the problem of evaluating the flux ~h- ds. Let dr be the line segment joining the two 
points ( x,t ) and (x+dx,t+dt) (see Fig. 2.2(a)). Let 7^ = (n x ,n t ) be a unit normal to dr. Then 

dt _ dx 


n x = ± 


^(dxj 1 + (dt) 2 


n, = T 


V (dx) 2 + {dt) 1 


(2.7) 


where {dx) 2 + (dt) 2 * 0 is assumed. The upper and lower signs in Eq. (2.7) correspond to the two 
senses ofT^ Let ds be the surface element with the end points (x,t) and (x+dx,t+dt). Then 


ds 


^ (dx) 2 + (dr) 2 7? = ±(dt , -dx) 


Eqs. (2.4) and (2.8) imply that 


■7) — > 

h ■ ds = ± 


( a u-V.^)dt-iidx 


= ±t-dt 


( 2 . 8 ) 


(2.9) 


where 


~t d - (~u , au-[L^) , dr ^ (dx , dt) (2.10) 

It may be shown that the upper (lower) signs in Eqs. (2.7) - (2.9) should be chosen if the 90° 
rotation from 7? to dr is in the counterclockwise (clockwise) direction. 

Let T be a simple closed curve in £ 2 - Let (x,t) and (x+dx,t+dt) be two points on T. Let 7? be 
the outward normal to T at the point (x,t) (see Fig. 2.2(b)). Then the upper (lower) signs in Eqs. 
(2.7) - (2.9) should be chosen if dr points in the counterclockwise (clockwise) direction of T. Let 
AT be a segment of r. Then Eq. (2.9) implies that 

f “T^ f C C ' > “7^ 

J „ h ds = l£ dr (2.11) 

J Ar ■'Ar 

where the notation c.c. indicates that the line integration should be carried out in the 
counterclockwise direction. 


Let OPQRS be the parallelogram depicted in Fig. 2.1(b). Let 

J(PQ) d - the flux of leaving OPQRS through the line segment PQ. (2.12) 

Similarly, we define J(QR), J(RS) and J(SP). Let T be the boundary of OPQRS. Then Eqs. (2.6), 
(2.10) and (2.1 1) may be used to obtain 


where 


Ar 

J(PQ) = ~ 


>m = f 


Ar 

J(RS) = -y- 


J(SP) = 


Ax 


d+t)y 1 


0-t)7; 


-(l+x)Y " 


-(1 — t) y7 


+ (1 -t 2 -8)ja" 

- (1 - x 2 -5)-y a" 

+ (l-x 2 + 5 )^a1 
4 J 

- ( l - x 2 + 8 ) -y a" 


(2.13) 

(2.14) 

(2.15) 

(2.16) 


and 


def ( a — b ) Ar 
Ax 


(2.17) 


8 V so 

(Ai) 2 

Two comments may be made about Eqs. (2.13) - (2.16): 
a. These equations are consistent with the local conservation relation 

J(PQ) + J(QR) + J(RS) + J(SP) = 0 
and 


(2.18) 


(2.19) 


b. The influence of parameters a and b on the fluxes leaving OPQRS through its four edges is 
expressed through a single parameter x. As a result, the use of the moving mesh depicted 
in Fig. 2.1(a) does not increase the complexity of the expressions on the right sides of Eqs. 
(2.13) - (2.16). The meaning of x is a subject to be discussed in Section 3. 

To proceed, let 

/i ( 0 ) C/» d - y-J(PQ) , / 2 (O) 0» d - - 7 - J(QR) ( 2 . 20 ) 
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( 2 . 21 ) 


^ — f-J(RS) , / 2 (/) 0',«) d - -~KSP) 

Ax Ax 


In other words, f[ 0) (J,n) and fj 0) (J,n), respectively, are the normalized fluxes leaving CE(J,n) 
through its "future right" and "future left" edges. Similarly, fi n (j,n) and / 2 (/) (/,«)» respectively, 
are the normalized fluxes entering CE(J,n ) through its "past left” and "past right" edges. For 
simplicity, a normalized flux will be referred to simply as a flux. Thus, the first two normalized 
fluxes may be referred to as the outgoing fluxes while the last two normalized fluxes as the 
incoming fluxes. These two pairs of fluxes form two column matrices, i.e.. 




./ 2 (O) 0» . 




'f} n (j,n) ' 

.fPij.n) . 


( 2 . 22 ) 


Also we define 


and 


<7i0» ^ y] 


<720'.«) ^ -^-a" 




def 


<7i0» 

?20'.«) 


With the aid of Eqs. (2.20) - (2.24), Eqs. (2.13) - (2.16) can be rewritten as 

7 i0) (j,n) = A {0) 7t(j,n) 
and 

/ (/) 0'.«) = A (/) 4 > (/,«) 

where A^° ^ and are the matrices defined by 


A (0) dg 

1 + x 1 - x 2 - 5 


1 -x -(1-x 2 -5) _ 

A (/) d g 

1 +X -(1 -x 2 + 8) 


1-T 1 -x 2 + 8 


Through out this paper, it will be assumed that 

1 - x 2 + 8 * 0 

Thus , the inverse of exists. We have 


(2.23) 

(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 
(2.29) 
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[A<'>r‘ = { 


1 

1 -X 


1 

1 +T 


1-t 2 + 5 1 - x 2 + 8 -I 


From Eq. (2.26), one obtains 

?(/.») = [AV>r l T n (J,n) 

Substituting Eq. (2.31) into Eq. (2.25), one has 

y {0) (j Un) = nf (I) (j,n) 

where is the matrix defined by 

a ^ A (0) [A (/ >r ! 


(2.30) 


(2.31) 

(2.32) 

(2.33) 


Let (o tm , /, m = 1, 2, be the elements of the matrix Q. Then Eqs. (2.27), (2.30) and (2.33) imply 
that 


(l-x 2 ) + 8 „ (l+x)(l-x 2 ) 

1 - x 2 + 8 ’ Wl2 - 1 - x 2 + 8 

_ ( 1 — x) ( 1 - x 2 ) „ _ -x(l-x 2 ) + 5 

1 - x 2 + 5 ’ 0)22 " 1 — x 2 + 8 



A result of Eqs. (2.34) and (2.35) is 

2 

2>/m = 1 , m = 1,2 

/ = i 

Since Eq. (2.32) is equivalent to 

// 0) (i.«) = E . i = U2 

m - 1 


Eq. (2.36) may be used to prove that 


/i (O) 0» + /2 (O) 0.«) = /i (,) 0» + A n (j.n) 


(2.34) 

(2.35) 

(2.36) 

(2.37) 

(2.38) 


i.e., the sum of the outgoing fluxes is equal to the sum of the incoming fluxes. From Eqs. (2.20) 
and (2.21), it is easy to see that Eq. (2.38) is equivalent to Eq. (2.19). 

With the above preliminaries, the current marching procedure may now be defined by using 
the interface flux conservation relation. Explicitly, this relation requires that (see Fig. 2.3) 

f{ n (j+'A,n+'/i) = /! (O) 0» . fP{j+ l /i,n+Vi) = / 2 (0) (/+!.«) (2.39) 


and 
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/i (/) 0>+l) = ti 0) (j- l A,n+'A) , / 2 (/) 0>+l) = fi° ] (j+ l / 2 ,n+V 2 ) (2.40) 

where j = 0, ±1, ±2, • • • ; n =0, 1, 2, • • ■ . Because of the above relations, a single arrow is 
drawn across an interface (see Fig. 2.3) to represent both the flux entering and the flux leaving 
this interface. Let n > 0 be a fixed integer. Let the outgoing fluxes of the conservation 

elements at the time level n be given. According to Eq. (2.39), the incoming fluxes 
fi n (j+ i A,n+ 1 /i) and (j +'A, n +‘A) of CE (J+lA.n+lA), respectively, are equal to the outgoing 

flux fl 0) (j,n ) of CE (j,n) and the outgoing flux fj 0) (j+l,n ) of CE(/+l,n). Thus all the 
incoming fluxes of the conservation elements at time level n+ l A are known. Since Eq. (2.37) 
remains valid if the indices j and n, respectively, are replaced by j+ l A and n+ l A, one has 

2 

fi (0) (J+V 2 ,n+'A) = 2 ®imfm ) <J+ l 6,n+ x A) , /= 1, 2 (2.41) 

m- 1 

The outgoing fluxes of the conservation elements at time level n+ l A can be evaluated in terms of 
the known incoming fluxes by using Eq. (2.41). Similarly, with the aid of Eq. (2.37) and (2.40), 
the incoming and outgoing fluxes of the conservation elements at time level n+ 1 can also be 
evaluated. This procedure can continue for time levels n+ 3/2, n+2, ■ • ■ . The following 
comments are made to provide more details about this procedure: 

a. The outgoing fluxes ft (0) (J, 0) at time level n = 0 may be evaluated by using Eq. (2.25) if 
the coefficients 0) at time level n = 0 are given. 

b. The coefficients co, m are considered as given constants in the marching procedure. Let j 
and n be a pair of fixed integers. Then the outgoing fluxes // 0 ) (/' +V 2 , n +'A), 1= 1,2, may 
be evaluated in terms of the incoming fluxes f^ r) (J+V 2 ,n+ l A), m = 1,2, by using Eq. 
(2.41). This evaluation requires four multiplications and two additions. However, the 
operation count can be reduced if the following alternative procedure is adopted. Let 
fl° ) (j+ l A,n+ 1 A) be evaluated by using Eq. (2.41). This requires two multiplications and 
one addition. fj°^(J+ l A,n+ l A) is then evaluated by using the conservation relation 

fi 0) (j+ l A,n+V 2 ) - f} l) (j+'A,n+'A ) + f^Hj+V^n+'A) - fi 0) (J+Vi,n+'A) (2.42) 

which can be obtained from Eq. (2.38) by replacing j and n, respectively, with j+ l A and 
n+'A. The last evaluation requires only one addition and one subtraction. Thus evaluation 
of ff 0 ^(j+V 2 ,n+'A), 1= 1, 2, requires only two multiplications, two additions, and one 
subtraction. Thus the operation count of the current scheme is five for each j per half time 
step. 

c. The principal variables involved in the marching procedure described above are incoming 

and outgoing fluxes. However, for any pair of given integers or half integers j and n, 
if (/,/t) can always be evaluated in terms by using Eq. (2.31). 
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More relations connecting to the marching procedure may be derived. To proceed, note that 
Eqs. (2.39) and (2.40), respectively, can be written as 

7 (/) 0 '+%,»+V4) = I + 7 (O) 0',«) + L/ O> 0+1,«) (2.43) 

and 

T n (j,n+\) = l+7 i0) (j ! - l An+ l A) + L? ( ° ) (j+'A,n+'/ l ) (2.44) 

where 


' 1 

0 ' 

, and 

L * 

' 0 

0 ' 

0 

0 



0 

1 


I + ¥ 


are projection matrices [p.l 16, 5], For any pair of numbers Cj and cj, we have 
I + 


(2.45) 


' c 1 ' 


’ Cl " 


~C 1 


' 0 


= 


, I_ 




. C 2 . 


. 0 j 


. C 2 . 


. C 2 . 


(2.46) 

Substituting Eq. (2.43) into the equation obtained from Eq. (2.31) by replacing j and n, 
respectively, with j+Vi and n+Vi, one obtains 

7fo'+ l /2,«+^) = [a (/) ] _i [i + 7 (0) o» + i-? (0) o+i.«>] 

With the aid of Eq. (2.25), Eq. (2.47) implies that 

t(j+'^n+Vi) = Q + ^(/,") + Q-^O'+l.") 

where 

Q+ % [A (/) r‘ I + A (0) , Q_ % [A^]" 1 L A<°> 

Multiplying Eq. (2.47) from the left by A (0 } and using Eqs. (2.25) and (2.33), one has 

/ ^ (j+V 2 ,n+V?) = n+y^ 0 ^(j,n) + £2_7* 0) (/'+l,/i) 

where 

def T 0 

^ qi + = 

©21 o 

Similarly, with the aid of Eq. (2.44), one can obtain 

~$(j,n+\) = Q + ~cf (j -Vi,n+Vi) + Q_ ^(J+Vi,n+'A) 
and 


Q_ ^ QL = 


0 ©12 
0 ©22 


(2.47) 


(2.48) 


(2.49) 


(2.50) 


(2.51) 


(2.52) 
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(2.53) 


T 0) (j\n+\) = n + y { °\j-V2,n+V2) + n_T 0) (j+V 2 ,n+'A) 


As a result of Eqs. (2.48) and (2.52), one has 

!t(J,n+\) = [ Q+ ] 2 ~<t(j—l,n) + [Q + Q_ + Q_Q + + [Q_ (2.54) 

Similarly, Eqs. (2.50) and (2.53) may be used to obtain 

/ O) 0'."+ 1) = [O+?T 0) (j-hn) 

+ [n + ^ + n_Q + ]y i0) (j,n) + [O-] 2 T 0) (j+l.n) (2.55) 


In both Eqs. (2.54) and (2.55), a column matrix at time level n+1 is expressed directly in terms of 
three column matrices at time level n. Note that, with the aid of Eqs. (2.23), (2.27), (2.30), (2.45) 
and (2.49), Eq. (2.54) can be explicitly expressed as 


yj 



X + 


5(1-t) 
1 - x 2 + 5 


(l+x)y*_i + 


(1 -t 2 -5) 


Ax a"_j 


1 -x 
-x 2 + 5 


(1-t 2 )y?- 


x(l — x 2 — 5) 


Ax a" 



5(1 +x) 
1 -x 2 + 5 


( 1 — X ) Y/+I 


(1 -x 2 -5) 

4 


Ax a" + ] 


(2.56) 


and 


Axa? +1 = - 4 - 


x + 


5(1-t) 
1 - x 2 + 5 


1 — x 2 
1 - x 2 + 8 


Y/-i + ( 1 - x) 


1 - x 2 - 5 
1 -x 2 + 5 


Ax a"_i 


1 -X 2 

2 

1 -x 2 + 8 

- 


4xyy + (1 -x 2 -5)Axa" 


x - 


5(1 +x) 
1 - x 2 + 5 




(2.57) 


At this juncture, it is noted that the marching procedure described earlier is constructed by 
using Eqs. (2.50) and (2.53). One may construct alternative procedures by using Eqs. (2.48) and 
(2.52), or Eq. (2.54), or Eq. (2.55). However, these alternatives are less efficient than the original 
procedure. This is because (i) the matrices 0 + and fl_, respectively, have only two surviving 
elements, and (ii) the original procedure can take advantage of the conservation relation Eq. 
(2.38). Since the current numerical method is developed on the principle of flux conservation, it 
is natural that the most efficient marching procedure is the one that uses incoming and outgoing 
fluxes as the marching variables. 
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To further clarify the differences between the current scheme and other explicit schemes, this 
section is concluded with the following remarks: 


a. The marching step in the Lax-Wendroff scheme in which is updated, i.e., Eq. (1.7), is 
derived with the assumption that u is a constant along a straight line with dx/dt = a. This 
assumption generally breaks down if u satisfies Eq. (2.2) instead of Eq. (1.3). If the 
diffusion term in Eq. (2.2) is comparable to the convection term, the error caused by this 
assumption may be substantial. Nevertheless, the marching step Eq. (1.7) or one of its 
variants is used in many generalized Lax-Wendroff schemes that are used to solve Eq. 
(2.2). This marching step generally is followed by another in which «" +1 is obtained by 
using the conservation relation (see Fig. 1.2(a)) 


,n + l 


Ax - u" Ax + 


„ \* + ' A 
au j+ l A ^ V 0^ )j 


)j+'A 


At 


aWj - |x( 


du \ft+'A 
dx 


\ n-r /7. 

)j-'A 


At = 0 


(2.58) 


where and respectively, are the 

dx J dx J 


finite-difference approximations of 


du/dx at the mesh points (j+ l /i,n+Vi) and (j-'/i,n+Vi). Generally, these approximations 
may be expressed in terms of the mesh values of u at the time level n+Vi. 


b. The assumption that u is a constant along a straight line with dx/dt= a may be avoided if 

one lets uJ+'a . 7 = 0 , ±1, ±2, • ■ ■ , be lagged behind by one half time step, i.e., = 

u'j+'A’j = 0, ±1, ±2, • • • . Here may be obtained by interpolating the given values of 
«y. j = 0, ±1, ±2, • ■ • . Obviously, an explicit conservative scheme may be formed by 
combining this new assumption with Eq. (2.58). 

c. The errors caused by the assumptions mentioned in (a) and (b) generally are considered as 
the penalty one pays for using an explicit conservative time-accurate scheme. One may 
avoid this penalty by using either an implicit scheme or an explicit nonconservative time- 
accurate scheme (e.g., the MacCormack scheme). The current scheme is an exception to 
the above common wisdom. It is explicit, conservative and time-accurate, yet constructed 
without relying on the assumptions mentioned in (a) and (b). 

d. In Section 1, the pure -convection discrete conservation relation Eq. (1.9) was cast into an 
integral form (i.e., Eq. (1.14)) similar to the conservation law Eq. (1.1). For \i * 0, one may 
be tempted to repeat the same feat by again assuming Eq. (1.12) but replacing Eq. (1.13) 
with 
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(2.59) 


d - ( au(x,t ) - ,u(x,t )) 

OX 

However, because du(x,t)/dx = 0, Eq, (1,14) again is equivalent to Eq. (1.9). It will not be 
equivalent to a convection-diffusion conservation relation in the form of Eq. (2.58). 

A desire to cast the discrete conservation relation into an integral form is a strong 
motive behind the current development. As a matter of fact, the integral form Eq. (1.18) is 
one of the basic building blocks of the current marching scheme. It is our belief that a 
conservative scheme that can be cast into an integral form not only is easier to interpret but 
also provides a more realistic simulation of the conservation laws. 

e. In the current scheme, u{x,t) is approximated by u(x,t) which is defined in Eq. (2.6). For 
(x,t) e CE "(J,n), n(x,t) is determined by two independent parameters y * and a" which, 
respectively, represent « and du/dx at the point (x",t n ). The extra parameter a" accorded 
to the current scheme allows a more precise specification of the discrete initial conditions. 
It also provides more leeway for &(x,t) to simulate a rapidly varying function u(x,t), as 
often occurs across a shock or within a boundary layer. 

f. According to Eqs. (2.32) and (2.33), the determination of/ (O) 0',«) in terms off U) (j,n) 
requires the inversion of the matrix A^ r \ As a result, the current scheme is locally implicit. 
As will be shown in Section 4, the appearance of the factor (1-t 2 + 5) in the 
denominators of two elements in [A (/) ] _1 (see Eq. (2.30)) has a positive effect on stability. 

g. Let 

1 — T 2 — 5 * 0 (2.60) 

Then [A (0 ^] -1 and exist. Let 

0 + ^ ^ Lfi' 1 (2.61) 

By using Eq. (2.32) and the interface flux conservation conditions, we obtain the time- 
reversal counterpart of Eq. (2.55), i.e., 

7 {0) (J,n) = [6_] 2 / O) 0-M+ 1) 

+ [£2 + O_ + Q_Q + ]7 (O) 0'»«+l) + [O^/’Wn+D (2.62) 

Eq. (2.62) states that the discrete variables associated with a conservation element at time 
level n can be determined by those of its three closest conservation elements at time level 
n+ 1. Contrarily, in a typical classical scheme, e.g., the Lax-Wendroff scheme Eq. (1.11), a 
discrete variable at time level n may depend on all the discrete variables at time level n+ 1. 
Note there may be several solutions ofT^O'.rt) for a given set of 
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O'+I.h+I), if the current scheme is generalized to solve a nonlinear 


PDE, e.g., the Burgers’ equation [p. 154, 3]. 
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3. THE DYNAMIC SPACE-TIME MESH 


The main purpose of this section is to explore the concept of a dynamic space-time mesh and 
the need for a unified treatment of physical variables and mesh parameters. Specifically, we will 
demonstrate that stability and accuracy of a numerical calculation may be improved if the space- 
time mesh is allowed to evolve with the physical variables such that the local convective motion 
of physical variables relative to the moving mesh is kept to a minimum. To simplify the 
discussions, again we consider only Eq. (1.3) or Eq. (2.2). Also the coefficients a and ft, and the 
mesh parameters b. At, and Ax are assumed to be frozen at their local values. 

The parameter x defined in Section 2 plays a central role in the following discussions. As a 
result, its role as the Courant number for a moving mesh will be established immediately. 

In Fig. 2.1(a), Q and S, respectively, denote the mesh points (j,n+Vi ) and (j,n- l A). The point 
T is on time level n-V 2 with TQ being in the direction of convection. Hereafter, by definition, a 
line segment in space-time is said to be in the direction of convection if dx /dt = a along this line 
segment. Note that the direction of convection is identical to the characteristic direction of Eq. 
(2.2) only if ft = 0. Points T and S are on the same time level and separated by a spatial distance 
( a—b )At. The parameter x is the ratio between this distance and Ax. In the case where b = 0, 
i.e., the moving mesh is reduced to a stationary mesh, the spatial distance between T and S is 
reduced to aAt. As a result, x is reduced to the ordinary Courant number v. For this reason, the 
parameter x may be considered as the Courant number for a moving mesh. 

To further explore the significance of x, again we consider Eq. (1.3) and the Lax-Wendroff 
scheme which solves it. If the moving mesh depicted in Fig. 2.1(a) is used, then this scheme may 
be expressed as 


and 


_ I 
J + ' A ~ 2 


( 1 + X ) M* + (1— X) Uj + i 


«; +1 = uj-% 


..n +V6 


.. n -f Vil 

- «>-V4 J 


(3.1) 


(3.2) 


As in the derivation of Eqs. (1.7) and (1.10), Eq. (3.1) is obtained through the use of backward 
characteristic projection and linear interpolation while Eq. (3.2) represents a flux conservation 
relation over the parallelogram PUVR shown in Fig. 2.1(a). When b = 0, x = v and Eqs. (3.1) and 
(3.2), respectively, are reduced to Eqs. (1.7) and (1.10). For this reason, the original Lax- 
Wendroff scheme may be viewed as a special case of the scheme defined by Eqs. (3.1) and (3.2). 
As will be shown later, a moving mesh relative to a coordinate system may become a stationary 
mesh relative to another coordinate system. As a result, a scheme defined by Eqs. (3.1) and (3.2) 
with any value of b may be reinterpreted as the original Lax-Wendroff scheme if it is viewed 
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from another coordinate system. This will provide an alternative (and more satisfying) proof for 
Eqs. (3.1) and (3.2). 

Note that v is a function of Ax, At, and a while x is a function of Ax, At, a, and b. The extra 
independent variable b of the function x corresponds to the extra degree of freedom introduced as 
a result of allowing the mesh to be "moving" relative to the coordinate system. It will be shown 
immediately that the time-step size limitation associated with the original Lax-Wendroff scheme 
may be removed by taking advantage of this added freedom. 

According to the von Neumann analysis, the amplification factor of the scheme defined by 
Eqs. (3.1) and (3.2) is 

A L W (Q) = 1 - t 2 (1 -cos0) - /xsinG (3.3) 

where 0 is the phase angle variation in Ax of a plane-wave component. Eq. (3.3) implies that the 
stability condition is |x| <1, i.e., 

4 ,£ T^T (a * b) <3 - 4 > 

Let Ax be held constant. Then Eq. (3.4) implies that the stability bound for At becomes greater as 
| a—b | becomes smaller. Since a and b, respectively, are the convection velocity of the physical 
variable u and the velocity of the moving mesh, a -b is the velocity at which u is convected 
relative to the moving mesh. In this paper, a — b and \a—b\, respectively, may simply be 
referred to as the relative convection velocity and the relative convection speed. As a result, one 
may say that the time step size limitation associated with the Lax-Wendroff scheme is due to the 
existence of a nonzero relative convection speed. 

A large relative convection speed and thus a severe time-step size limitation, may result from 
an indiscriminate use of a stationary mesh. For the current case in which a is a constant, this 
limitation may be eliminated completely by using a moving mesh with b = a. Even if a is a 
function of u, x, and t, the above discussion suggests that the time-step size limitation may be 
reduced sharply if the mesh is designed such that the local relative convection speed is kept to a 
minimum. 

At this juncture, we introduce another interpretation for the parameter x, i.e., it is the product 
of the relative convection velocity ( a-b ) and the mesh aspect ratio (At/Ax). If only the 
stationary mesh (i.e., b= 0) is allowed, then x can be made smaller only by making At /Ax 
smaller, a move very costly in computational effort. On the other hand, if a moving mesh is 
allowed, x may be made smaller by reducing | a — b | . 

Next we study the dependency of accuracy on the relative convection speed \a -b | . Note 
that A l w (0) = 1 when x = 0. Thus the numerical dissipation and dispersion vanish when x = 0 
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[pp. 93-94, 3]. Moreover, Eq. (3.2) implies that the numerical solution does not vary along a j 
mesh line when x = 0. Since a j mesh line is also a characteristic line of Eq. (1.3) when x = 0, the 
numerical solution coincides with the analytical solution at the mesh points. Since x — » 0 as 
(a-b) —> 0 when At/ Ax is held constant, the above observations suggest that the accuracy of the 
scheme defined by Eqs. (3.1) and (3.2) may also be improved by reducing \a—b\. 

Since A LW ($) = when x = ±l, the dissipation and dispersion also vanish when x = ±l. 
Again Eqs. (3.1) and (3.2) imply that the numerical solution is exact when x = ±l (Note: The 
diagonal RU (PV) of the parallelogram PUVR depicted in Fig. 2.1(a) is in the characteristic 
direction of Eq. (1.3) if x = 1 (x = -l)). Thus the numerical solution of Eqs. (3.1) and (3.2) will 
be highly accurate if one uses a mesh with |x| <1 and |x| is very close to 1 everywhere. 
However, since | x | = 1 is on the verge of instability, this strategy of obtaining accurate solutions 
may not be practical when the coefficient a is not a constant. 

In order to further explore the significance of x and ( a-b ), in the following, we will study 
the transformation properties of several equations and parameters under a Galilean 
transformation. This study will also provide a systematic way to obtain the form of a classical 
finite-difference scheme over a moving mesh. 

To proceed, we consider the Galilean transformation: 


x - x-b*t and t' = t ( 3 . 5 ) 

where b* is any real constant. Physically, (x',t') represents a coordinate system moving with the 
velocity b * relative to the coordinate system (x,t). Assuming that the mesh is fixed in space- 
time, Eqs. (2.3) and (3.5) imply that the coordinates x' and t' for the mesh point (j,n) are given 
by 

x' = x'j jAx + n b'At and t' = t ,n d =f nAt (3.6) 

where 


b' 


def 


b-b * 


(3.7) 


is the velocity of the moving mesh relative to the new coordinate system. In this paper, a 
parameter defined with respect to the new coordinate system is denoted by a prime. Immediately, 
we have 


Ax' ^ x '% j - x ,n j = Ax 


At' ^ t' n+1 


- t' n = At 


Also, Eq. (2.2) may be rewritten as 


(3.8) 



(3.9) 


du' 

dt' 


+ a’ 


du' 

dx 


- M- 


av 

dx' 2 


= 0 


where 

a' ^ a- b* , \Jl' \l (3.10) 

and u' is a function of x' and t' such that 

u'(x',t') = u(x,t ) (3.11) 

An immediate result of Eqs. (3.7) and (3.10) is 

a'-b'-a-b (3.12) 

Thus the relative convection velocity (a -b) is invariant under the Galilean transformation Eq. 
(3.5). Also Eqs. (3.8), (3.10) and (3.12) imply that 

t' V = [a-b)t* m t 5- M ijx^. . _ J 

Ax' Ax (Ax') 2 (Ax) 2 

Moreover, it may be shown that, for any (x,t) e CE "(J,n) 

u(x,t) = f % a'Jix'-x']) + P'"(t' -f'") + Yj (3.14) 

where 

a'; & a] , p'7 % p; + b'a* , Yj ^ Yj (3.15) 

With the aid of Eqs. (3.10) and (3.15), one concludes that a'a'" + p'"=0 if and only if 
a a" + P" = 0. 


From Eqs. (3.8), (3.13) and (3.15), one concludes that all the parameters and variables that 
appear in Eqs. (2.56) and (2.57) are invariant under the Galilean transformation Eq. (3.5). This 
property will be used to simplify the discussion given in Section 6. 

Let b* =b. Then b' =0, and thus Eqs. (3.6), (3.8) and (3.13) imply that 

x'j=jAx' , t' n = nAt' (b*=b) (3.16) 

and 


x = x' 


a'Ar' 

Ax' 


(**=*) 


(3.17) 


From Eqs. (3.16) and (3.17), one concludes that (i) the mesh is stationary relative to the primed 
coordinate system, and (ii) x simply becomes an ordinary Courant number when it is expressed in 
terms of the primed parameters. Moreover, as a result of (i), a classical finite-difference scheme 
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may now be expressed relative to this new coordinate system in its traditional form. As an 
example, the L/D-F scheme fp. 161, 3] for solving Eq. (3.9) may be expressed as 


f n+\ 


- U 


2 At' 


'"-1 i/ ,n j - h " 1 „'»+! _ 

, “ i+i ~ “ j-i _ / a j±l zJL±±ZJLi u l 

2Ax' ( Ax ' ) 2 


+ a 


(b*=b) 


= 0 


(3.18) 


Since the mesh is fixed in space-time, Eq. (3.1 1) implies that u'j = u* for any (J,n). With the aid 
of Eqs. (3.8) and (3.10), Eq. (3.18) may be rewritten as 




- U 

2 At 


n-\ 


+ ( a-b ) 


Uj+l ~ 1 

2Ax 




«/+i 


+ ul 


-t 

(Ax) 2 




- u 


n-1 


= 0 (3.19) 


This is the form of the L/D-F scheme when the mesh and coordinate system used are those 
depicted in Fig. 2.1(a). As a result, when a stationary mesh (b = 0) is replaced by a moving mesh 
(b * 0) without changing the coordinate system, the only modification required in the form of the 
L/D-F scheme is to replace the coefficient a with (a-b). This is also true for other classical 
schemes solving Eq. (1.3) or Eq. (2.2). Note that the Courant number v should be replaced by the 
parameter x as the coefficient a is replaced by (a-b). This is consistent with the fact that Eqs. 
(1.7) and (1.10), respectively, are converted into Eqs. (3.1) and (3.2) when v is replaced by x. 

In conclusion, the previous discussions suggest that a reduction in the relative convection 
speed | a -b \ may improve stability and accuracy, and reduce dissipation and dispersion of 
numerical calculations. Since the appearance of wiggles near a discontinuity is a result of 
numerical dispersion, the wiggles may also be reduced by reducing the relative convection speed. 
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4. STABILITY, DISSIPATION AND DISPERSION 

4.1 Preliminaries 

In this section, the current numerical scheme will be studied using a discrete Fourier analysis. 
Specifically, we assume the initial periodic conditions: 0) =~tf(j+K, 0), ;' = 0, ±1, ±2, 

where K is an integer > 3. With the aid of Eq. (2.54), by induction, it may be shown thatT^ (/',«) 
is periodic at any time level n, i.e., 

= ~i?(j+K,n) 0=0, ±1,±2, ••• , n = 0, 1, 2, ••• ) (4.1) 

With the aid of Eqs. (2.54) and (4.1),^(/» can be expressed explicitly as a matrix function of j, 
n, K and the initial-value matrices if(/, 0), l = 1, 2, 3, • • • , K- 1. The stability, dissipation and 
dispersion of are then studied by using this functional relation. According to Eqs. (2.25), 

(2.26) and (2.48), the other matrix variables, including !?(j+V 2 ,n+ l A), 7^ (/.«). y^(j,n), 
(J+ l A,n+ l / 2 ) and \j + ] /i,n+Vi), may be considered as functions of ~c?(j,n). Their 

behaviors, therefore, may be inferred directly from those of !?(J,ri). 

Since the current Fourier analysis also serves as the basis of an error analysis to be presented 
in Section 5, the following development will include materials that are needed there. 


(4.2) 

(4.3) 


Z* (*,*' = 0, 1, 2, • • ■ , K - 1 ) 

r=o 

where 5**' is the Kronecker delta symbol. As a result, 

K - 1 

?0» = Z ?(*.«) <t>f U =0, ±1, ±2, ■ • • , n =0, 1, 2, • ■ ■ ) 


£(*.«) d ~ zV(/.«)t>f* ) (* = 0, 1,2, -,K-\ , n=0, 1,2, ••• ) 

/=o 


(4.4) 


(4.5) 


(4.6) 


To proceed, let 


d =f -^Lr exp [2nijk/K] 


/-C T 

, *=0, 1,2, • • • , AT— 1 ) 


0 =0, ±1, ±2, ■ 

(j)j A) are periodic and orthonormal, i.e., 

<1 >f = <t>5+AT O =0, ±1, ±2, • ■ • , k = 0, 1 , 2, ■ • • , K-\ ) 
and 
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Hereafter, unless specified otherwise, it is assumed that it = 0,1,2, ■ • • , AT-1; 
j = 0, ±1, ±2, • • • ; and n= 0, 1, 2, 

Furthermore, let 

Q d ef | 2tc k/K if K/2>k >0 

1 2 n(k-K)/K if K-\>k> K/2 
and 

Q (6) & e~ im Q + + e im Q. ( jc > 0 > - * ) 

Note that 9*, £ = 0, 1 , 2, • • • , AT-1 , are deliberately defined such that 

Jt > 0* > -Jt 

Also, unless specified otherwise, hereafter we assume that Jt > 0 > -jc. Substituting Eq. (4.5) into 
Eq. (2.54), and using Eqs. (4.3), (4.4), (4.7) and (4.8), one has 

t( k ’ n+ 1) = [<2(0*)] 2 ^f(*,/») (4.10) 

i.e., the amplification matrix for any k is the square of the matrix (2(0*). Combining Eqs. (4.2), 
(4.5) - (4.7) and (4.10), it may be shown that 

± £ IQ Ok)] 2 ” z' e i{J ~ mk ?(/. 0) (4.1 1) 

A k = 0 /=0 

i.e., the matrices q*(j,n) are determined uniquely if the initial-value matrices ^(/,0), 
/ = 0, 1 , 2, ■ • • , K-l , are given. 


(4.7) 

(4.8) 

(4.9) 


With the aid of Eqs. (2.27), (2.30), (2.45), (2.49) and (4.8), it can be shown that 
cos(9/2) - i Tsin(0/2) -i ( 1 -x 2 -5)sin(0/2) 

1 -x 2 -5 


( 2 ( 0 ) = 


/ ( 1 -x 2 )sin(0/2) 
1 -t 2 + 5 


(4.12) 


1 -t 2 + 8 


[ cos(0/2) + i x sin(0/2) ] 


Let 


11(0) d - 6 cos(~) - /x(l-T 2 )sin(-|-) 

Then the eigenvalues of Q (0) are 

CT (0) dg r|(9) ± ^[r|(9)] 2 + ( 1 -t 2 ) 2 - S 2 
± 1 - x 2 + 8 

In this paper, the principal square root is defined such that 


(4.13) 


(4.14) 
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— S the phase angle of its polar form > — — 

By applying the von Neumann analysis to Eq. (2.54), it may be shown that the amplification 
factors are the eigenvalues of [ Q (0)] 2 , i.e., [a + (0)] 2 and [ct_(0)] 2 . 

To proceed further, note that, for each 0, the matrix Q (0) is either nondefective or defective 
[p.353, 6], If Q (0) is nondefective, the Jordan form (2(0) of £2(0) and its powers [(2(0)] 2 , 
[Q(0)] 3 . • • • may be chosen as [p.362, 6] 

[ [o + (0)] m 0 1 

[Q(0)] m = m = 1, 2, 3, • • • (4.15) 

L 0 [a_(Q)r J 

On the other hand, if Q (0) is defective, we have [p.362, 6] 

a + (0) = a_(0) 
and 

[o + (0)] m mto+tf)]'"- 1 ' 

0 [a_(0)] m . 

According to Jordan’s theorem [p.362, 6], for each 0 (tc>0>— k), there exists a nonsingular 
matrix G (0) such that 

(2(0) = G(0) Q(0) [G(0)] _1 (4.17) 

Note that matrix G (0) is not unique. It can be shown that a matrix G(0) can also convert Q (0) to 
Q(Q) if and only if G(0) = G (0)^(0) where ¥(0) is (i) an arbitrary 2x2 nonsingular diagonal 
matrix if Q(0) has two distinct eigenvalues, or (ii) an arbitray 2x2 nonsingular matrix if Q(Q) is a 
multiple of the identity matrix, or (Hi) a 2x2 nonsingular upper-triangular matrix with identical 
diagonal elements if £2(0) is defective. The above comments are useful in a later discussion. 

Substituting Eq. (4.17) into Eq. (4.1 1), one arrives at 

?(/.«) = V e m G(0* ) [£2(e*)] 2 "4 (4.18) 

k = 0 

where the column matrices are defined by 

g j 

4 ^ ^ [G(B k )]~ l £^ 0 ‘?(/,O) (4.19) 

K 1=0 

Note that there is a one-to-one relation between the set of column matrices c*. k = 0, 1, 2, • • • , 
K- 1, and the set of column matrices 'if (/, 0), j = 0, 1, 2, • • • , K- 1. As a matter of fact, one has 


[Qm m = 
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(4.20) 


Let 


and 


t(j.0)= Z e ij6k G(Q k )-? k 

k = 0 


G(Q) = 


gum 

g 21 0) 


^ 12(B) ■ 
# 22 ( 0 ) 


(4.21) 




Sn(8) ' 
g 21 (0) 


?-(0) 


g 12 ( 0 ) 
. £22(0) 


(4.22) 


With the aid of Eqs. (4.15), (4.16), (4.21) and (4.22), Eq. (4.17), which is equivalent to Q (0)G(0) 
= G(0)G(0). implies that (i)~g + (Q) is an eigenvector of (2(0) with the eigenvalue o + (0), (ii)t-(Q) 
is an eigenvector of Q (6) with the eigenvalue a_(8) if (2(0) is nondefective, and (Hi) j2(0)^_(0) 
= £+(0) + a_(0)£>_(0) if Q(0) is defective. 


Let Q (0*). k = 0, 1,2, • • • , K- 1, be nondefective or defective with a + (0*) = a_(0*) = 0 
(Note: [Q (0*)] 2 = 0 in the latter case). Then Eq. (4. 1 8) is reduced to 

K-\ 

[a + (0*)] 2n c* + g > + (0*) + [o_(0*)] 2n c*_^_(0*) } (4.23) 

*=o 


where c k+ and c*_, respectively, are the upper and lower elements of the column matrix c*. 
Several comments may be made relating to Eq. (4.23): 

a. The influence of the initial-value matrices ^(/, 0) on t(J> n ) is expressed through the 
coefficients c k+ and c*_. 

b. Let 

t ± (/>.*) e‘ j0t [o ± (0*)] 2n c k± t±(B k ) (4.24) 

Then, for each k,~tf(j,n) = t+(J,n,k) or l?(j,n) = l?_(J,n,k) is a particular solution of Eqs. 
(2.54) and (4.1). The general solution given in Eq. (4.23) is the sum of these particular 
solutions. 


c. With the aid of Eqs. (2.45), (4.19), (4.21) and (4.22), and the fact that c k+ and c*_, 
respectively, are the upper and lower elements ofc*, it can be shown that 

Ctt = G (0*) I ± % = ± G (0*) I ± [G (0*)]- 1 *£ e~ U6t f(l, 0) (4.25) 

A /=o 

Let the eigenvalues of (2(0*) be distinct. Then, as noted earlier, matrix G(0*) which 
converts Q (0*) into (2(0*) (in this case, (2(0*) is a diagonal matrix) can be replaced by 
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Q(9k) = G(0 i ) K F(0 i ) where ¥(0*) is an arbitrary 2x2 nonsingular diagonal matrix. Since 
I+, I_, ^(0*) and ['F(0*)] _1 are diagonal and thus commute among themselves, 

G(0*) I ± [G(Q k )T l = G(0t) m) I± t'I'O*)]- 1 [G(9 k )r l = G(0*) I ± [G(9 k )]~ l (4.26) 

Combining Eqs, (4.25) and (4.26), one concludes that the matrices <^±^±(6*) are invariants 
under the transformation G (0*) -> G(0*) if the eigenvalues of Q (0*) are distinct. 

To interpret the particular solutions defined in Eq. (4.24), we introduce the functions P ± (0) 
such that 

[o±(0)] 2 = | o ± (0) | 2 e ‘ P±(0) and 7t> p ± (0) > -jc (4.27) 

Note that p + (0) ( p_(0) ) is uniquely defined by Eq. (4.27) if a + (0) * 0 ( o_(0) * 0 ). Also we 
define 


d J$ u Ax 

a±(0) ^ b- 


0 


def 1° I <J±(0) I 

M±0) ^ rrr 


At 

2 


/ 0 \ 2 a 


(0*0) 

(4.28) 

(0*0) 

(4.29) 


and 


p±{x,t,9) J 




2— i B+(0) — 

|o±(0)| * e - A ' 


if 0*0 
if 0 = 0 


By using Eqs. (2.3) and (4.27) - (4.30), it may be shown that 

e ij6t [a ± (0*)p = P±(x],t H ,9 k ) 


(4.30) 


(4.31) 


According to Eqs. (4.24) and (4.31), ~?±(j,n,k) is the product of p±(xj,t n ,9 k ) and c*±g±(0*)- 
Since the latter is independent of j and n, the behaviors of 7?±(j,n,k) are governed by the former. 


For any 0 such that re > 0 > — tc and 9*0, u =p±(x,t,B) is a plane-wave solution of the 
convection-diffusion equation 


3 u .... 3 u ... 

IF + IFF - fe (6) 


3 2 u 

dx 2 


= 0 


Also, the wavelength of this solution is given by 

m v 


(4.32) 


- 33 - 



As a result of the above observations, for each k * 0 (i.e., 0* * 0), the particular solution 
~ct(j,n)=lt±(j,n,k) may be referred to as the plane-wave solution with the numerical covection 
speed a±(9*), the numerical viscosity ^(0*) and the wavelength X(0*). Also since 0 O =0 and 
p±(x,t, 0) is independent of x, one may say that the particular solution =~ct±(j,n, 0) has an 
infinitely long wavelength. 

In this paper, the marching procedure defined by Eqs. (2.54) and (4.1) is said to be stable if 
and only if, for any integer K >3 and any specification of the matrices c*. k = 0, 1, 2, • • • , K-\ 
(i.e., any specification of the initial-value matrices ^(/, 0). See Eqs. (4.19) and (4.20)), the 
elements of the column matrices j = 0, ±1, ±2, • • • , remain bounded as n —> +« with 

the parameters x and 5 being held constant (i.e., At and Ax being held constant — if one assumes 
that a, b and p are constants). The readers are reminded that the term "stability" referred to in 
Lax’s equivalence theorem has a meaning different from what we define here (see Section 6). 

Because G(Q k ), k= 0, 1, 2, • ■ ■ , K- 1, are nonsingular, Eq. (4.18) implies that the marching 
procedure is stable if and only if, for any integer K 1 3, the elements of the matrices [2(0*)] 2 ' 1 , 
k = 0, 1, 2, ■ • • , AT-1, remain bounded as n — » +°° with the parameters x and 5 being held 
constant. According to Eqs. (4.15) and (4.16), this implies that stability occurs if and only if, for 
any K > 3 and any k =0, 1, 2, ■ ■ ■ , K- 1, 

max { |a + (0*)| , |a_(9*)| } < 1 if Q (0*) is nondefective (4.33) 

and 


|a + (0*)| < 1 if Q (0*) is defective (4.34) 

Our study of dissipation and dispersion will be limited to the case in which each matrix Q (0 t ) 
is either nondefective or defective with o + (0*) = o_(0*) = 0. From Eqs. (4.23), (4.24), and (4.27) 
- (4.32), one concludes that (i) for any k * 0, the dissipation of ~<t ± (j,n,k) may be measured by 
J4±(©jfc), (ii) the dispersion of the general solution ~<?(j,n ) may be determined by the 
distribution of a±(0*), k- 1,2, • • • , K- 1. 

As a final note of this subsection, we will point out a remarkable similarity between the forms 
of a ± (0) and the amplification factors A ± (0) (see Eq. (B.10)) of the L/D-F scheme. Let 
1 -x 2 > 0. Then both the numerator and denominator on the right side of Eq. (4.14) may be 
divided by 1 — x 2 . Asa result, we have 


o±(0) = 


5cos(-|) - /xsin(y) ± 


V 


[6cos(-|) - t xsin(y)] 2 + 1 - 5 2 


1 + 5 


(4.35) 


where 5 8/(1 - x 2 ). A comparison between Eqs. (4.35) and (B.10) reveals that the expression 
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on the right side of Eq. (B.10) may be converted to that on the right side of Eq. (4.35) if 5/2, t 
and 0, respectively, are replaced by 5, x and 0/2. In making this comparison, the reader should 
keep in mind that the amplification factors for the current scheme are [o+(0)] 2 and [o_(0)] 2 , 
rather than o + (0) and o_(0). Note that, if 1 - x 2 < 0, the sign "±" on the right side of Eq. (4.35) 
should be replaced by "T". 

This completes the preliminaries. A discussion of two special cases, i.e., ( i ) 5 = 0 and (ii) 
x = 0, will precede the investigation into the general case in which both 5 and x may not vanish. 

4.2 The Special Case With 5 = 0 

Eqs. (2.29) and (4.14) coupled with the assumption 6 = 0 imply that (i) x 2 * 1, and (ii) 

o ± (0) = - / x sin(0/2) ± -L 1 ~ T I Vi -x 2 sin 2 (0/2) (4.36) 

1 — x z 

In the following, Eq. (4.36) will be used to study (a) stability and dissipation, and (b) dispersion. 

(a) Stability and Dissipation : 

We have 

1 if |xsin(0/2)| < 1 

M e >l = j . |,_ t 2, (4.37) 

| -xsin(0/2)± -L— r L V x 2 si n 2 (0/2) - 1 | if |xsin(0/2)| > 1 

In the case where x 2 > 1, there exist a K and a k (K-\ >k>0) such that |xsin(0*/2) I > 1. Thus 

max { |a + (0*)| , |o_(0*)| ) > |xsin(0*/2)| > 1 (4.38) 

Combining Eqs. (4.33), (4.34), and (4.38), one concludes that the current marching procedure is 
not stable if x 2 > 1 and 5 = 0. 

1^1 x 2 < 1. The Eq. (4.36) implies that a + (0) * o_(0) for any 0. As a result, the matrices 
Q(Q), k > 0 > — jc, are nondefective. According to Eq. (4.37), we also have |a ± (0)| = 1, n > 0 > 
-k. It follows from Eq. (4.33) that the marching procedure is stable if x 2 < 1 and 5 = 0. 
Moreover, since |o±(0*)| =1, k = 0, 1, 2, • • ■ , K- 1, the particular solutions defined by Eq. 
(4.24) will not dissipate as n increases. Thus the numerics reflects faithfully the physics of pure 
convection. This contrasts sharply with the Lax-Wendroff scheme Eq. (1.11) which is 
numerically diffusive even though it is also a conservative scheme. 

(b) Dispersion : 


- 35 - 



Let x 2 < 1. Then Eq. (4.36) implies that 

[o ± (0)] 2 = g 2i Sin -1 [ t sin(8/ 2) ] (439) 

Since Go = 0 and [a ± (0)] 2 = 1, the particular solutions defined in Eq. (4.24) are independent of j 
and n if k = 0. As a matter of fact, the term with k = 0 on the right side of Eq. (4.23) is reduced to 
a constant column matrix co«- g*+(0) + Co_^_(0). Thus this term is ignored in the following 
discussion of the dispersion of the general solution defined in Eq. (4.23). 

Since the range of Sin -1 is (-71/2,71/ 2], one concludes that 

f > Sin _1 [xsin(0/2)] > -y if t 2 < 1 (4.40) 

As a result, a comparison between Eqs. (4.27) and (4.39) reveals that 

K > p ± (0) = T 2 Sin -1 [ x sin(0/ 2) ] > -n (4.41) 

Substituting Eq. (4.41) into Eq. (4.28) and using Eq. (2.17), one concludes that 

(e*0) (4.42) 

Thus 


2±(0) = a if x = 0 and 0*0 (4.43) 

i.e., when x = 0, all the particular solutions which appear on the right side of Eq. (4.23) except the 
one with k = 0 are "convected" with the same velocity a. In this case, dispersion is completely 
absent. 

In the case where 1 > x 2 > 0, Eq. (4.42) may be expressed as 


where 


fl±(0) = a-x[lTF(x,0)] 


Ax 

At 


(l>x 2 >0; 0*0) 


ev t m def Sin 1 (xsin(0/2)) 
1 ’ } x(0/2) 


(1 > x 2 > 0 ; 0*0) 


It is an exercise in calculus to show that 

1 > F(x,0) > 2 /t 1 

By using Eqs. (4.44) and (4.46), one obtains that 


(4.44) 


(4.45) 


(4.46) 
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(1 >x 2 >0; 0*0) 


(4.47) 


0 < [a-fl + (0)]/( x^) < 1-- 

' At K 

and 

2 Ar 

1 + — < [a- a_(0) ] / (x — ) <2 (l>x 2 >0; 0*0) (4.48) 

Eqs. (4.47) and (4.48) state that, on the real line, the distance between any a + (0) and the physical 

2 

convection velocity a is less than (1 ) |x | Ax/ At while the distance between any a_(0) and a is 

71 ~ 

2 

greater than (1 + — )|x|Ax/Ar and less than 2 1 x | Ax/ At. Thus the dispersion, measured by the 

maximum spread between a and any a+(0) or a_(0), is less than 2 1 x | Ax /At. As | x | decreases, so 
does the dispersion. Recall that the same conclusion was also reached in Section 3 for the Lax- 
Wendroff scheme. Moreover, since the maximum of the spread between a and a+(0) is less than 
the minimum of the spread between a and any <z_(0), a particular solution defined by taking the 
upper sign in Eq. (4.24) will be "convected" at a velocity closer to the physical convection 
velocity a than a particular solution defined by taking the lower sign in Eq. (4.24). For this 
reason a + (0) and ^+(0), respectively, may be referred to as the principal eigenvalue and 
eigenvector of matrix £2(0) while o_(0) and g*_(0), respectively, the spurious eigenvalue and 
eigenvector of £2(0). This designation may be extended to the case 0 = 0 even though a±(0) are 
undefined at 0 = 0. Similarly, a particular solution ^?±(j,n,k) will be referred to as a 

principal (spurious) solution if the upper (lower) sign is chosen. 


4.3 The Special Case With x = 0 


For this special case, the physical variable u has no convective motion relative to the moving 
mesh. Also Eq. (4. 14) is reduced to 


_ 5 cos(0/2) ± V i _ 5 2 sin 2 (0/2) 

a ±(°) = T71 


(4.49) 


Eq. (4.49), coupled with (i) 5 > 0 and (ii) cos(0/2) > 0 if jc > 0 > -k, implies that 
5cos(0/2) ± Vl - 5 2 sin 2 (0/2) 1 2 


l ° ± ( 0)| 2 =\ 


1+5 


5-1 

5+1 


< 1 


< 1 if |5sin(0/2)| < 1 


if 1 5 sin(0/2) | > 1 


(4.50) 


Moreover, with the aid of Eq. (4.27) and the additional definitions that (i) P + (0) 0 if a + (0) = 0 

and (ii) P_(0) 0 if a_(0) = 0, one concludes that 
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0 


if |5sin(0/2)| <1 


P±(9) = i 


± 2 Sin" 1 


V 5 2 sin 2 (0/2)~ 1 
5 2 - 1 


if |8sin(0/2)| > 1 


(4.51) 


To study the stability, note that a + (0) = o_(0) is a necessary condition for the defectiveness of 
matrix Q(Q), and it occurs if and only if 1 8 sin(0/2) | = 1. Asa result, it follows from Eqs. (4.33), 
(4.34), and (4.50) that the current scheme is unconditionally stable if x = 0. 


The dissipation and dispersion of the particular solutions defined by Eq. (4.24) generally may 
be studied explicitly by using Eqs. (4.28), (4.29), (4.50), and (4.51). This study is greatly 
simplified if one considers only the case in which 1 ^SsiO. For this special case, 
1 5sin(0/2) [ < 1 for any 0. According to Eq. (4.51), p±(0) = 0. n > 0 > -tc, i.e., the dispersion is 
absent. Moreover, by studying the extrema of the first expression on the right side of Eq. (4.50), 
it may be shown that 


1 > M 0)| 2 £ |- p | > | a _( 0)| 2 > ( 1 > 8 > 0 ) ( 4 . 52 ) 

According to Eqs. (4.24) and (4.52), the rate of dissipation per time step of any spurious solution 
is greater than or equal to that of any principal solution if 1 > 8 > 0 and x = 0. 


4.4 The General Case 

Assuming 8 > 0 and 1 — x 2 + 8 ^ 0, it is shown in Appendix A that the current scheme is 
stable if and only if t 2 < 1. This stability condition has the remarkable property that it is 
independent of p except that p > 0 is assumed. 

Assuming 8 > 0, it is shown in Appendix B that the stability region of the L/D-F scheme on 
the 8-x plane is the region defined by x 2 < 1, minus the two points (0,1) and (0,-1). This 
stability region is exactly identical to that of the current scheme. 

On the other hand, the stability conditions of all other classical schemes known to the authors 
are dependent on p. As an example, the stability region of the MacCormack scheme (see 
Appendix D) is depicted in Fig. 4.1. Obviously, the stability region of the MacCormack scheme 
is smaller than that of the L/D-F and the current scheme. The significance of this difference was 
discussed in Section 1. It will be further studied in sections 5 and 7. 

The dissipation and dispersion properties of the current scheme may also be studied for the 
general case in which both 8 and x may not vanish. This requires the use of Eqs. (4.27) - (4.29) 
and (A. 10). 
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5. ERROR ANALYSIS 


In this section, an error analysis technique is developed using the discrete Fourier analysis 
formulated in Section 4. Assuming smooth initial data, this technique enables us to predict, 
analyze and compare the numerical errors of the L/D-F, the MacCormack, and the current 
schemes for calculations involving hundreds or thousands of time steps. As will be shown, the 
results of this error analysis provide us with a theoretical basis for improving the accuracy of the 
current scheme. They will also be used to interpret the numerical results to be presented in 
Section 7. 

As a preliminary, the error analysis will be preceded by a discussion on a notable feature of 
the current scheme, i.e., the requirement to specify two sets of initial data involving the values of 
Yy and a°. 

Among the classical schemes, the initiation of the L/D-F scheme also requires the input of 
two sets of initial data, i.e., u® and u ) . Since only u® are given, generally uj are evaluated in 
terms of u® by using a starting condition, e.g., Eq. (B.l). Since the starting scheme is constructed 
with the aid of an one-sided difference approximation of a time derivative, it is one order less 
accurate than the main scheme. As a result, the accuracy of the L/D-F scheme may not attain the 
level that one would expect if only the main scheme is considered. 

In the current scheme, the initial data Yy and ay, respectively, will be identified with u® and 
(du/dx) j. In the case where u(x, 0) is smooth and known for all x on the initial line, both u® and 
(du/dx)® may be evaluated and used as the initial data for the numerical calculation. Generally, 
the extra set of initial data ( du/dx )® will allow a more accurate approximation of u (x, 0) and thus 
gives the current method an edge in obtaining more accurate numerical solutions. 

In the current error analysis, the accuracy of the MacCormack, the L/D-F and the current 
methods will be studied and compared assuming that only the initial data «y are given. For the 
current method, this means that (du/dx)® must be evaluated in terms of u®. Since du/dx at any 
point on the initial line t = 0 is the result of the differentiation of u along the initial line, (du/dx)® 
may be expressed in terms of u® without using a one-sided difference approximation. Thus, at 
least in principle, the accuracy of the current scheme may not be reduced as a result of the 
complication associated with the extra initial data (du/dx)®. This contention will be verified by 
the results of the following analysis. 

To proceed, let the initial data u®,j = 0, 1, 2, • ■ ■ , K - 1 be given. Any u®,j = 0, ±1, ±2, 
may then be determined by using the periodic condition u® = uy+* , j = 0, ±1, ±2, ■ • • . In view of 
Eq. (2.6), it is natural to assume that 

Yy = u® , j = 0, ±1, ±2, (5.1) 
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In order to determine ( du/dx )° in terms of w°, and also provide the initial values for a 
corresponding analytical problem (i.e., u(x, 0 ) for this analytical problem will be determined in 
terms of u]), a smooth periodic fimction I(x) will be formed by linearly combining K periodic 
exponential functions (see Eq. (5.5)) such that 

(«) I (/Ax) = ufj , y = 0, 1, 2 , ■ • • , AT -1 ( 5 . 2 ) 

and 

(/>) I(x+KAx ) = I (x) ( 53 ) 

As a result of Eq. (5.2), and the fact that ay is the spatial derivative in CE"(j, 0), we will assume 
that 


«/ = /(/Ax) , j = 0, ±1, ±2, • • • (5.4) 

where l(x) is the derivative of I (x) with respect to x. 

Given yj and o tj , the discrete solution to Eqs. (2.54) and (4.1) may be determined. Assuming 
m(x, 0) — /(x), the analytical solution to Eq. (2.2) with the periodic condition 
u(x+KAx,t) = u(x,t) may also be determined. The accuracy of the discrete solution may then be 
assessed by comparing it with the analytical solution. 

To construct I (x), note that the exponential functions 

//(x) ^ e' WCATA*) f / = 0 , ±i, ± 2 , • • • ( 5 . 5 ) 

form a basis for the function space of the functions that have period KAx and are of bounded 
variation over [O^TAx] [p.478, 2], Let the integer K* > 1 be defined by 

, def f (AT— 1 )/2 if A' is odd 

K ~ I (5.6) 

l K/2 if AT is even 

Then it may be shown that the wavelengths of the K functions 7/(x), / = K*- K + 1, K*-K + 2, 

• • ■ , K*, are longer than or equal to those of the other functions defined in Eq. (5.5). We assume 
that I (x) is a linear combination of these K functions. With the aid of Eqs. (4.7) and (5.2), it may 
be shown that 

AT — 1 j j 

/(*) = Z b k e ** 1 (5.7) 

k-0 

where 
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(5.8) 


b k 


& i S 1 e-' 70 * u? 
K 1=0 


k = 0, 1,2, ■■■ ,K- 1 


At this point, it should be emphasized that the function I(x) defined by Eqs. (5.7) and (5.8) 
generally is complex even if the given initial data u° are real. As a result, generally a° may be 
complex. This should not cause alarm since the discrete equation to be solved, i.e., Eq. (2.54), is 
a system of algebraic equations with real constant coefficients. For these equations, both real and 
imaginary parts of a complex solution are themselves solutions. As a result, in case that physics 
so dictates, only the real parts of the initial data and solution may be considered as physically 
relevant. Obviously, the above comments are also applicable to a differential equation with real 
constant coefficients like Eq. (2.2). 


To obtain the solution to Eqs. (2.54) and (4.1), note that a result of Eqs. (2.23), (2.24), (4.2), 
(4.4), (4.7), (5.1), (5.4), (5.7), and (5.8) is 

1 


1 ^ -t78, 


K 


£ <T"V(',0) = b k 


1 = 0 


iOk 

4 


(5.9) 


Combining Eqs. (4.25) and (5.9), one obtains 

c*±g*±(9*) = b k Jt±(Q k ) , k= 0, 1,2, ■ • • , K— 1 


(5.10) 


where 


#±<e) ^ G(0)i ± tG(0)r 1 


1 

10 


l 4 J 


By definition, (1+ + 1_) is the 2x2 identity matrix. Thus Eq. (5.1 1) implies that 


tf + (0) + #_(0) = f} a (Q) ^ 


1 

10 


L 4 J 


(5.11) 


(5.12) 


In this section, we assume that every (2(0*), k = 0, 1, 2, •••, K- 1, has two distinct 
eigenvalues. As a result, Eqs. (4.23) and (4.26) are applicable. With the aid of Eq. (5.10), Eqs. 
(4.23) and (4.24) imply that 

tV-n) = *Zt(j '.".*) (5-13) 

*=o 


where 
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~tf(j,n,k) ~tf + (j,n,k) +l?_(j,n,k) 

= b k e ijQk { [ a + (0*) J 2 " 7? + (0*) + [ a_(0*) l 2 " j}_(Q k ) } (5.14) 

According to Eq. (5.14 ),~i?(j,n,k) is composed of the principal and the spurious parts. At n = 0, 
the amplitudes of these two parts are b k Jf + (O k ) and b k j}_(Q k ), respectively. 

Let u (x, 0) = / (a:). Then the analytical solution to Eq. (2.2) with u (x+KAx,t) = u (x,t) is 


u(x,t) = u a (x,t) ^ X b k e 

k = 0 


K \ ( \ 2 • ^ 


(5.15) 


Ua(xJ,t n ) 


taij.n) & r , . 

Ax du a (x,t) 


(5.16) 


z = x" ,1 = t’ 


In view of Eqs. (1.17a), (1.17b), (2.23), and (2.24), ~<? a (j,n) may be considered as the analytical 
counterpart ofrf(j,n). Let 


A " {rt+ T' )0 

A a (0) - e 


(5.17) 


Z(J,n,k) & b k e ii9i [AMT%^ k ) 

Then Eqs. (5.12), (5.15), (5.16), (2.3), and (2.18) may be used to show that 

K - 1 

^0» = E Z(j>n,k) 

*=o 

Combining Eqs. (5.12), (5.14), and (5.18), one obtains that 

&tf(j,n,k) ~$(j,n,k)-lt a (J,n,k) = b k e ,j9k [i? + (/i,0*) + £?_(«, 0*)] 


where 


£±(n, 0) ^ { [ O±(0) ] 2n - [ A a (0) ]" ) 7}±(Q) 


(5.18) 


(5.19) 


(5.20) 


(5.21) 


Thusl?(j,n,k) =~<? a (j,n,k) if both [a+(0)] 2 and [o_(0)] 2 are replaced by A a (Q). For this reason, 
A„(0) may be referred to as the analytical amplification factor. Because ~(? a (j,n) is the analytical 
counterpart o H?(j,n), the numerical error of tf(j,n) may be measured by 
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(5.22) 


AT- 1 


A ^0» - ?(/.«) -^(/.n) = z 

k= 0 

The last equality sign follows from Eqs. (5.13), (5.19) and (5.20). 

In the following, will be studied assuming 

1 > x 2 and 5 > 0 

To proceed, we define 

©(0) d - cos(0/2) - i x sin(0/2) - a_(0) , ji> 0>-7 c 


(5.23) 


(5.24) 


In Appendix C, it is shown that (i) co(0) * 0 if n > 1 0 1 , and (ii) <d(jc) * 0 if either l -% 2 - 5 * o or 
0 > x > -1. In the following discussion, we assume that <y(0) * 0. Let 

t def ( 1 -x 2 )sin( 0/2) - def ( 1 -x 2 -5)sin(0/2) 

5,< ’ (i-x 2 +5)to(e) ’ fe<6) " 5© (5 ' 2S) 

Then it is shown in Appendix C that, for any 0 such that jc > 0 > -n and o + (0) * a_(0), one has 


(a) 


(b) 


and 


(c) 


l+5.(0)5 2 (0) * 0 

0^2(0) r 1 


^+(0) = 


l + 


l+5l (0)52(0) 


1 m . t[ ^' m 

!+5i (0)52(0) 


i'5i(0) 
i 5 2 (0) 


(5.26) 


(5.27) 


(5.28) 


Note that, for the special cases in which 5 = 0 or x = 0, (0) and £ 2 (0) are given by 

' sin(0/2) 


5i(0) = 


cos(0/2) + Vi -x 2 sin 2 (0/2) 

sin(0/2) 

cos(0/2) + V] _6 2 sin 2 (0/2) 


if 5 = 0 


if x = 0 


(5.29) 


and 
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( 1 -x 2 )sin(0/2) 
cos(0/2) + Vi -x 2 sin 2 (0/2) 


if 5 = 0 


C2O) = 


( 1 — 5 2 ) sin(0/ 2) 
cos(0/2) + Vi - 5 2 sin 2 (0/2) 


if T = 0 


(5.30) 


In view of the different roles played by the parameters 5 and x, the structural similarity among the 
expressions on the right sides of Eqs. (5.29) and (5.30) is indeed remarkable. 


According to Eq. (5.14) and the comment following Eq. (4.32), the wavelength of any 
particular solution n,k) is inversely proportional to |0* | if k *0 and Ax is held 

constant. This is also true for its analytical counterpart. Thus a particular solution with a smaller 
1 0*| is a slower-varying function of the index j. In the following, we will study &tf(J,n,k) 
assuming 1 0* | is small. Note that a general solution which is a slow-varying function of the 
index j generally is dominated by the particular solutions with small 1 0* | , i.e., the coefficients b k 
are very small except for those k's with very small 1 0* | . 

As a preliminary to the following use of the Taylor’s expansion, note that, according to Eq. 
(4.14), the current assumption o + (0) * o_(0) is valid if and only if 

CO) d - [r|(0)] 2 + (l-T 2 ) 2 -5 2 *0 (5.31) 

By assumption 1 > x 2 . Thus C(0) = ( 1 -x 2 ) 2 > 0 and <o(0) = 2 ( 1 -x 2 ) /( 1 -x 2 + 5) > 0. It 
follows that that there is a neighborhood of 0 = 0 on the complex 0- plane in which both V^(Q) 
and to(0) are nonzero analytical functions of 0 (Note: In obtaining the Taylor’s expansions for the 
following study, the functions involved may be considered as the hmctions of a complex variable 
0 ). 


By using Eq. (5.28), it may be shown that 


#- 0 ) = 


(l-X 2 )2-5 2 
32(l-x 2 ) \ 

x80 2 


8 ( 1 — x 2 ) 192 


7 l|r 0 3 + [ ( 1 + sx 2 ) -t- 1 ] e 4 + o o 5 ) j 

1 [(1+3x 2 )+ 3§2 j (1 ~ 5 f ) ]0 3 +O(0 4 ) 


(5.32) 


(1-x 2 ) 2 


Hereafter a quantity is denoted by O (Q 1 ) if there exists a constant C > 0 such that the absolute 
value of this quantity < C \ 0 1 1 for all sufficiently small 1 0 1 . Eq. (5.32) is reduced to 



#-( 6) = 


/Do 

~ ^ ( 1 192 } 93 + <9 ^ 4 > 


if 5 = 0 


(5.33) 


or 


ft-(0) = 


( 1 - 5 2 ) ( 1 + 3 5 2 ) o4 
768 


/ ( 1 + 3 5 2 ) n3 . „ /rt 4 


192 


e 4 + <9(0 5 ) 
e 3 + o(9 4 ) 


if x = 0 


(5.34) 


Note that Taylor’s expression of //+(0) may be obtained from that of 7/_(0) by using Eq. (5.12). 
An inspection of these expansions reveals that the upper and lower elements in 7/_(0), 
respectively, are smaller than those in //+(0) by (i) three orders and one order of 0 if x 5 * 0 or (ii) 
four orders and two orders of 0 if either 5 = 0 or x = 0. Thus, at n = 0, the spurious part of 
is much smaller than the principal part if j 0* | « 1. 

In the case where 1 0| c 1, 77 ± (0) may be approximated by the leading terms in their Taylor’s 
expansions. In the following, we will search for the approximations of 

[a ± (0)] 2n -U a (0)r 

which are valid for small 0 and large n. Note that the approximations which are valid for only 
small n are of little value since n generally is quite large in a typical numerical calculation. 


To proceed, let 


MB) * 

A a (0) 


(5.35) 


By definition, |e + (0)| is a measure of error when A a (Q) is approximated by [a + (0)] 2 . It may be 
shown that 


£+( 9 ) = 


'*A,(t, 5) ,3.5, l-Sx 2 . 


24 


0 3 + 


192 ^ l-x 2 


A,(x,6)-4x 2 ]0 4 + O(0 5 ) 


where 


Aj(x,5) 1 — x 2 — 


35 2 

l-x 2 


(5.36) 


(5.37) 


is an often-encountered expression in the current paper. 
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(0 

and 


Since 


[o + (e)] 2 "-[/i a (e)r = u a (0)]» {[l+MG)]"-!} 


(5.38) 


00 [ 1 + e+(0) ]" = 1 + n e+(0) if n | e + (0) | c 1 (5.39) 

we arrive at the important conclusion that 

[ o + (0) l 2 " - [A a (Q) ) n i n [A a (Q) ] n £+(0) if n |e + (0) | < 1 (5.40) 

Here the sign "=" is used to signal the fact that the ratio between the expressions on both sides of 
this sign is nearly equal to 1. 

Note that the assumption n | e+(0) | c 1 may be justified if one is only interested in the case in 
which the numerical calculation is accurate up to the time step n, i .e.,7f(j,l,k) ±~(f a (j,l,k), l = 1, 
2, 3, According to Eqs. (5.20) and (5.21), generally this requires that [o + (G*)] 22 = 

[Ai(9*)] , > l = 1,2, • ■ • , n. It is shown in Appendix C that the last n equations are valid if and 
only if n \ £+(0*) | cl. 

Let 


£_( 0 ) ^ 


[ q -( 9)] 2 

A 2 (x,5)A a (0) 


(5.41) 


where 


By definition, | e^(0) | 
may be shown that 


A 2 (t,5) ^ 


1 -t 2 -5 


(5.42) 


1 -t 2 +5 

is a measure of error when A 2 (t, 5) A a (0) is approximated by [a_(0)] 2 . It 


e_(0) = 2i T0 + (5/2 - 2t 2 )0 2 + O(0 3 ) (5.43) 

Eqs. (5.41) and (5.43) imply that [o_(0)] 2 A 2 (x,5)A a (0) as 1 0 1 — > 0. 

Because (i ) A 2 (t,5) < 1 if and only if 5(1 -x 2 ) > 0, and (ii) e_(0) = O(0), Eq. (5.41) implies 
that 


I [ q -( 9)] 2 1 

1 A a ( 0) I 


|A 2 (x, 5)[ 1 +e_(0)]| < 1 


(5.44) 


if 6 > 0 and |0| is sufficiently small (note: 1 -x 2 > 0 is assumed in Eq. (5.23)). Assuming Eq. 
(5.44), one may conclude that | [o_(0)] 2b / [A a (0)]" | c 1 and thus 
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[o_(0)] 2 ' , -[^ a (0)r = -[^(0)r 


(5.45) 


if n is sufficiently large. 

As an example, let x = 0.5 and 5 = 0.1. Then A 2 (x,6) = 0.585. Let 0 be sufficiently small 
such that | e_(0) | <0.1. Then | [a_(0)] 2 /A a (0)| <0.65. Thus |[o_(0)] 2 "/[A a (0)]" | < 1.812 x 
10 -4 if n = 20. 

In the following discussions, the integers n and k are such that the approximations given in 
Eqs. (5.40) and (5.45) are valid if 0 = 0*. Let b k *0. Let r t (n,k), i = 1, 2, denote the ratio 
between the t'-th element in the column matrix A 7?(j,n,k) and that in the column matrix lf a (j,n,k) 
[Note: the meaning of rfak) will be examined later]. Then Eqs. (5.18), (5.20) and (5.21) imply 
that 


r, •(«,*) = r i+ (n,k) + ri_(n,k) , i = 1,2 


(5.46) 


with 


, , x def [^±( n >9*)li . 

r i± (n,k) =* =5 , i = 1, 2 

[A a (0*)] n [// a (0*)] ( - 


(5.47) 


Here [Z? ± («, 0*)],- and [7? a (0*)],, respectively, denote the z-th elements of the column matrices 
Z? ± (n,0*) and // a (0*). Hereafter, r i+ (n,k ) and r t _(n,k), respectively, will be referred to as the 
principal and spurious parts of r,(n,Jfc). By using Eqs. (5.12), (5.21), (5.32), (5.36), (5.40), and 
(5.45), one has 


r !+(«>£) = n- 


zxA^x.5) , 5 r i-5x 2 


24 


+ T^T [ 


192 1 1 -x 2 


A 1 (t,5)-4t 2 ]0* + O(0*) 


r\ ~(n,k) = - 


(l-x 2 ) 2 -5 2 J i x5 


32 ( l -x 2 ) 1 — x 




1 _2 \ . 35 2 (1-9x 2 ) 1q4 

( 1 -x 2 ) 2 


+ ^-[d + 3x 2 ) + 


r 2+ (n,k) = n- 


iz A,(x,5) 3 5 r 1 — 9 x 2 


24 


0* 4- 


192 l-x 2 


] Qi +0(Q 5 k ) | 

AKx.S) -4X 2 ] 0* + 0(0*) 


and 


r 2 -(n,k) = 


z'x 5 


2 ( 1 — x ) 


5T e * + 48 


( 1 + 3 x ) + 


.2 % . 35 2 (l-5x 2 ) 


d-T Z ) 


2 \2 


ei + 0(el) 


(5.48) 


(5.49) 

(5.50) 


(5.51) 


Note that, in the current paper, any term denoted by O (Q ! ) or O (0*) is independent of n. 
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Now we shall reexamine the meaning of r x (n,k) and r 2 (n,k). By using Eqs. (5.15), (5.17), 
(5.18), and (5.12), it may be shown that 

AT-l 

u a (Xj,t n ) = £ u a (j,n,k) (5.52) 

*=o 


where 


u a (J,n,k ) - b k e ij ° k [A a (d k )] n = the upper element in !? a (J,n,k) 
Let u (j,n,k) be the numerical counterpart of u a (J,n,k), i.e., 

u(J,n,k ) - the upper element in l?(j,n,k) 

Then it may be shown that 

u(j,n,k)-u a (j,n,k) 

= — ~ (b k * 0) 


u a (j,n,k) 


Similarly, it may be shown that 


«*,(*, V") & 


du a (x,t) 

dx 


X=X; y l = r 


K - 1 
*=0 


(5.53) 


(5.54) 


(5.55) 


(5.56) 


where 


Uax(j,n,k) ^ l -~b k e ij%t [A a (Q k )r 


Ax 


x 


the lower element 


in ~tf a (j,n,k) j 


The numerical counterpart of Uax(J,n,k ) is 


f N def 4 
Hx(j ,n,k) ^ — x 
Ax 


the lower element 


in ~?(j,n,k) j 


Obviously 


, HxO< n >k) - UaxQ.n.k) 

ri(.n,k) = — (b k * 0) 

Max (J,n,k) 


(5.57) 


(5.58) 


(5.59) 


The accuracy of the current scheme will be compared with those of the L/D-F scheme and the 
MacCormack scheme. In Appendices B and D, we study the numerical solutions of Eq. (2.2) 
generated by these two schemes. In these studies, again we assume periodic conditions, i.e., Eq. 
(B.2) and use the mesh depicted in Fig. 2.1(a). 



Let \i^{j,n,k) be the numerical counterpart of u a (j,n,k ) obtained by using the L/D-F scheme 
(seeEqs. (B.31) and (B.33)). Let 


, IV def Uj.(j’ n ’k) - u a (j,n,k) 
r L (n,k ) ^ — — {b k * 0) 

u a (j,n,k) 


(5.60) 


Then it is shown in Appendix B that 

r L {n,k ) = r L+ (n,k) + r t _(«,fc) 


with 

r L+ (n,k) = «| ^^-9* + -^-(1-t 2 )(1--|-5 2 )0* 

+ [9x 2 5 2 ( l-x 2 ) + 3x 4 5 + 4x 2 - 8 2 +2]0* + 0(0*) 

96 2 


(5.61) 


(5.62) 


and 

r L -(n,k) = |(1--j-)|x 2 0?--^(1-2x 2 )0^ 

— -j— [4x 2 (4— 9x 2 ) + 3 8 2 ( 15 x 4 — 12 x 2 + 1 ) ] 0 4 + (?(0|) 1 (5.63) 


r L+ (rt,k), and r L _(n,k), respectively, may be referred to as the principal and spurious parts of 
r L (n,k). 


Let u M (j,n,k) be the numerical counterpart of u a (j,n,k ) obtained by using the MacCormack 
scheme. Let 


,,, def H»(j^,k)-u a (j,n,k) 

r u(n,k) ?r ; (b k * 0 ) 

u a (j,n,k) 

Then it is shown in Appendix D that 

r M (n,k) = nl 0* + ± [5( i _ 6x 2 ) -6x 2 ( l-x 2 )]0^ + O(0 5 k ) 


(5.64) 


(5.65) 


The numerical errors of the current scheme, the L/D-F scheme, and the MacCormack scheme 
will be studied and compared using Eqs. (5.46), (5.48), (5.49), (5.61) - (5.63), and (5.65). Note 
that the parameter ri(n,k ) and the associated equations, i.e., Eqs. (5.50) and (5.51), have no 
counterparts in the last two classical schemes. 
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According to Eqs. (5.48) and (5.49), the principal and the spurious parts of r x (n,k) are of the 
same order of 8*. Because the principal part is linear in n while the spurious part is independent 
of n, generally, one expects that the principal part will become dominant as n increases. This 
conclusion also applies to r L (n,k). However, it may not apply to r 2 (n,k ) because its principal 
part is n O (0*) while its spurious part is O (0*). 

By comparing Eqs. (5.48), (5.49), (5.62), and (5.63), one concludes that the principal and 
spurious parts of r x (n,k), respectively, are one order of 0* smaller than those of r L (n,k). As 
shown in Appendix B, the difference reflects the fact that the current scheme is more accurate 
than the UD-F scheme in both the initial-value specification and the main marching procedure. 

According to Eqs. (5.48), (5.49), and (5.65), r M {n,k), and both the principal and the spurious 
parts of r x (n,k) are in the same order of 0*. Moreover, one may observe that: 

a. For any x * 0, the ratio between the leading term in the principal part of r x (n,k) and the 
leading term in r M (n,k) approaches 1/4 as 5 — » 0. 

b. The leading term in the spurious part of r j (n,k) generally will be small for small 5. 

The above observations coupled with the fact that the principal part would be dorminant for a 
large n lead us to conclude that the current scheme generally will be more accurate than the 
MacCormack scheme by a factor of 4 when 5 is small and the initial condition is smooth. 

In the above discussion, the general constraints on x and 5, respectively, are 1 -x 2 > 0 and 
5 > 0. For the MacCormack scheme, these constraints must be further tightened by stability 
consideration (see Fig. 4.1). In the following, we will discuss the additional constraints required 
to annihilate the leading term in the principal part of each of the parameters r x (n,k), r 2 (n,k), 
r L (n,k), and r M (n,k) (Note: r u {n,k ) has only the principal part). This annihilation obviously will 
lead to a sharp improvement in the accuracy of the schemes under consideration. 


Let x = 0. Then all the leading terms in the principal and spurious parts of r x (n,k), r 2 (n,k), 
r t (n»k), and r u (n,k) vanish. As a matter of fact, it may be shown that 


r \(n,k) = n 


r 2 (n,k) = n 


8 ( 1 - 3 5 2 ) a4 . „, n5 , (l-5 2 )(l + 35 2 )„4 


02 + 0(0^) 


5(l-35 2 ) a 4 . — 5 \ . (1+3 8 2 ) a2 


02 + 0(92) 


92 + 0(9*) +-' Vo ^ 92 + 0(92) 


i(n,k) = n 8(4 19 2 ^ } 92 + 0(9*) - +0&1) 


(5.66) 


(5.67) 


(5.68) 


and 
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r u (.n,k) = n 


(5.69) 


48 

Note that: (i) The second leading terms on the right sides of Eqs. (5.62) and (5.63) also vanish 
when t = 0, and (ii) when x = 0, the ratio of the leading terms in the remaining r x +(n,k), r L+ (n,k), 
and r M {n,k) approaches 1:4:4 as 5 — > 0. Again we emphasize that, for a stationary mesh (i.e., b = 
0), x = 0 occurs only when a = 0. However, for any a * 0, x may be annihilated if one uses a 
moving mesh with b = a. 

According to Eq. (5.65), the leading term in r u {n,k ) also vanishes when x 2 = 1. However, 
since x 2 = 1 occurs either outside or on the stability boundary of the MacCormack scheme (see 
Fig. 4.1), the strategy of improving the accuracy of the MacCormack method by choosing the 
parameters At, Ax and b such that x 2 = 1 may be impractical in reality. This is particularly true 
for the more general case in which the convection speed may be a function of the dependent 
variable u. 


According to Eqs. (5.48) and (5.50), the leading terms in ri + (n,k ) and r 2+ («,/!:) also vanish 
when A] (x,8) = 0. Combining A] (x,5) = 0 and the assumption 1 -x 2 > 0, one obtains the optimal 
condition 


1 - x 2 = V 3 5 

Combining Eqs. (5.46), (5.48) - (5.51) and (5.70), it may be shown that 

x 2 (l-x 2 ) 


(5.70) 


r\(n,k) = n 


48^3 


et+orn 


(1-X 2 ) 

48 


^ e * + u' (1 ~ 3x2)0 * +o(0 * ) 


(5.71) 


and 


r 2 (n,k) = n 


x 2 ( 1 — x 2 ) 
48^3 


0 * + 6 >( 0 *) 


1 x 


1 


2V3- 0 * + ^ (1_x2)0 * +O(0 * ) 


(5.72) 


Eq. (5.70) represents a parabola on the 8-x plane. Since the segment of this parabola with 
8 > 0 lies entirely within the stability boundary of the current method, no stability constraint may 
prevent us from improving the accuracy of the current method by choosing At, Ax, and b such 
that 1 -x 2 = VJ 5 . 


Because 8 > 0, Eq. (5.70) is equivalent to the parametric equations: 


- 51 - 



(5.73) 


x = tan(y/2) (-£>^>-5.) 
and 

5 = [1 -tan 2 (\j//2)] (y>V>-y) 

Eqs. (5.73) and (5.74) imply that 

def X "^3" , . K K 

Re - g- = ~y tan ^V) ( y > V > -y ) 

With the aid of Eqs. (2.17) and (2.18), one also has 

Re = 1 £-&) A* 

4p 

As a result, Re may be referred to as the mesh Reynolds number. 


(5.74) 


(5.75) 


(5.76) 


In the case where ( a - b ), p and Ax are given, Re and y may be determined using Eqs. (5.76) 
and (5.75). Subsequently, the values of x, 5, and At when A^x.8) = 0 may be determined, 
respectively, by Eqs. (5.73) and (5.74), and the relation Ar = (Ax) 2 8/(4p). For later reference, 
these particular values of x, 5 and At, respectively, will be denoted by x 0 , 5 q and (At) 0 . 

In the case where p, At, and Ax are given, 5 may be determined using Eq. (2.18). If 
5 > 1 /VJ, Eq. (5.70) has no solution. If 1 /^3 > 5 > 0, then Eq. (5.70) implies that A ( (x,5) = 0 if 

x = x± ±Vi _ V3 5 (5 .77) 

Subsequently, the values of ( a-b ) when A!(x,5) = 0 may be determined by using the relation 
(a-b) = x Ax/At. 


Eqs. (5.49), (5.51), and (5.63) were obtained assuming Eqs. (5.45) and (B.43). The last two 
equations, in tum, assume 5 > 0. In the following, we consider the special case in which 8 = 0. 

Eqs. (4.37) and (5.17) imply that 


As a result, 


|o_(0)| = | A a (0) | =1 if x 2 < 1 and 8 = 0 


[q-(0)r-[A a (0)]" 

[A„m n 


<2 if x 2 < 1 and 8 = 0 


(5.78) 


(5.79) 


With the aid of Eqs. (5.12), (5.33), (5.47), and (5.21), one concludes that 



|r i_(n,*)[ < 


(5 = 0) 


(5.80) 


( 1-t 2 )(1+3t 2 ) n4 


384 


©2 + 0(0*) 


and 


|r 2 _0a)| ^l (1+ 2 4 — 0l + O(Ql) (5 = 0) 


(5.81) 


Eqs. (5.48) and (5.50) are applicable even if 5 = 0. In that case, they, respectively, are reduced to 

/ X ( 1 — X 2 ) a3 


ru.(n,k) + n 


24 


0* + 0(0*) 


(5 = 0 ) 


(5.82) 


and 


r 2 +(n,k) = n 


£ X ( 1 - X 2 ) 
24 


el + ool) 


(5 = 0) 


(5.83) 


A comparison between Eqs. (5.80) and (5.82) reveals that the spurious part of r\(n,k ) is 
negligible compared with its principal part if 5 = 0. Thus 


r\ (n,k) = r \+(n,k) = n 


/ T ( 1 — T 2 ) 
24 


Ql+0(Q 5 k ) 


(5 = 0) (5.84) 


Since r 2 -(n,k ) < O (0 2 ) and r 2+ (,n,k ) = nO (0*), the relative weights of the principal and spurious 
parts of r 2 (n,k) are dependent on the relative magnitudes of 0* and n if 5 = 0. 


By using an argument similar to that leading to Eqs. (5.80) and (5.81), it is shown in 
Appendix B that 


\r L _(n,k)\ < — 0* - 


X 2 „ 2 X 2 (4-9x 2 ) „ 4 ■ s\ ,r\5 


24 


0 * + 0 ( 0 *) (5 = 0 ) 


if x 2 < 1 and 5 = 0. Eq. (5.62) is applicable even if 5 = 0. In that case, it is reduced to 

r2- 


r L +(n,k) = n 


M±-L±Ql + 0(Q 5 k ) 


(5 = 0) 


(5.85) 


(5.86) 


Eqs. (5.85) and (5.86) imply that the relative weights of the principal and spurious parts of r L (n,k) 
are dependent on the relative magnitudes of 0* and n if 5 = 0. 

By comparing Eqs. (5.84), (5.86), and the reduced form of Eq. (5.65) for the special case 
5 = 0, i.e., 
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r M (n,k ) = n 


(5 = 0) (5.87) 


IT( 1-T 2 ) 3 

6 9 *~ 


t 2 (1-t 2 ) 


e^+oof) 


One concludes that the ratios of the leading terms in r 1+ (n,*), r L+ {n,k), and r M (n,k) are exactly 
1:4:4 if 8 = 0. 


Because comparison of the accuracy is meaningless without considering the operation counts 
of the schemes being compared, we conclude this section with a comparison of the operation 
counts of the current and the MacCormack schemes. It is shown in Section 2 and Appendix D 
that, for each j, (i) it requires four multiplications, four additions, and two subtractions for the 
current scheme to advance one full time step, and (ii) it requires five multiplications and four 
additions for the MacCormack scheme to advance one time step. 



6. CONSISTENCY AND THE TRUNCATION ERROR 

In Section 5, we studied the question of how well a discrete solution to Eq. (2.54) and (4.1) 
approximates its analytical counterpart. In this section, we will investigate the circumstances 
. under which an analytical solution may "satisfy" Eq. (2.54). This investigation amounts to a 
study of the consistency and the truncation error of the current scheme. This study will assume a 
uniform mesh like that depicted in Fig. 2.1(a). The analysis in this section will be further 
simplified by assuming b = 0, i.e., the mesh is stationary with respect to some coordinate system 
(x,t). According to a discussion given in Section 3, the last assumption may be made without any 
loss of generality. Note that consistency and the truncation error are properties to be evaluated at 
each point within the computational domain. The above uniform-mesh assumption is tantamount 
to freezing the mesh parameters at their local values. Also, the assumption b = 0 is tantamount to 
introducing a local coordinate system (x,t) such that the mesh is stationary with respect to this 
system at the local point under consideration. 

Before proceeding, we will discuss a general limitation on the ability of an explicit scheme to 
solve a convection-diffusion problem accurately. As an example, consider Eq. (2.2) (in this 
section, unless specified otherwise, p > 0 in Eq. (2.2)) over a domain with d>x> 0 and t > 0 (see 
Fig. 6.1). Let the initial data u(x, 0) (d>x>0), and the boundary data u(0,t) and u(d,t) 
(f > 0) be given. Let u(P 0 ) and u (P 0 ), respectively, denote the values of analytical and discrete 
solutions at a fixed point P 0 . Since a characteristic of Eq. (2.2) is represented by / = constant, the 
domain of dependence of u(P 0 ) is the union of AB, BC, and CD. In other words, u(P 0 ) is 
dependent on all the initial data, and the boundary data with t < t 0 . Assuming that the discrete 
solution is generated by an explicit solver, then the domain of dependence of u(P 0 )> contrarily, 
will include only a subset of the mesh points located on AB, BC, and CD. As an example, 
consider the MacCormack scheme (see Eq. (D.3)). If the mesh point (j,n) is not on or immediate 
next to the boundary, then uj +l is determined by u"_ 2 , j, uj, u"+i and u" +2 . As a result, the 
domain of dependence of a (P o) includes only the mesh points on EB, BC, and CF. Because (i) 
the mesh points that lie on AB but not EB and those that lie on CD but not CF do not belong to 
the domain of dependence of u(P 0 ), and (ii) for a fixed point P 0 , the lengths of AE and FD are 
proportional to the ratio At /Ax, one may conclude that, as At, Ax -» 0, the discrete solution 
( considered as a function of At and Ax) can not converge to its analytical counterpart unless 
At/Ax -4 0. 

As another example, we consider Eq. (2.2) with +°° > x > and t > 0 (see Fig. 6.2). Let 
u(x, 0) = u 0 (x) where u 0 (x) is a given smooth function with period d, i.e., u 0 (x+d) = u 0 (x). The 
solution to the above problem also is smooth and has period d in the x-direction. The domain of 
dependence of u(P 0 ) is the entire initial line. Assuming that the discrete solution is generated by 
the MacCormack scheme, then the domain of dependence of n(P 0 ), contrarily, includes only the 
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mesh points on EF. If to is small enough, the length of EF will be even less than the period d of 
the initial data. Because the length of EF is inversely proportional to the ratio At/Ax, one may 
conclude that, as At, Ax -> 0, the discrete solution can not converge to its analytical counterpart 
unless At/ Ax -> 0. Note that this conclusion was also mentioned by Fritz John on p.l 1 1 in [7], 

Contrarily, in the case of pure convection equation Eq. (1.3), the domain of dependence of an 
analytical solution at (x,t) is a single point (x-at , 0) on the initial line. As a result, At/Ax -> 0 
is not required for the convergence of an explicit-scheme discrete solution to its analytical 
counterpart. 

Because an explicit-scheme discrete solution to Eq. (2.2) will not converge to its analytical 
counterpart as At, Ax -> 0 without imposing the additional condition At /Ax 0, in general one 
would not expect that, in the limit of At, Ax —> 0, the nodal values of an analytical solution to Eq. 
(2.2) will satisfy the explicit scheme which is satisfied by the discrete solution at all At and Ax. 

In this paper, by definition, a discrete scheme is said to be "strongly consistent" with a PDE 
when the nodal values of any smooth solution of the PDE satisfy the discrete scheme (i.e., the 
truncation error = 0) in the limit of At, Ax — » 0, regardless how the mesh is refined (particularly 
how At/Ax behaves as At, Ax -» 0). According to the observation given in the last paragraph, it 
should be an exception rather than the norm for an explicit scheme to be strongly consistent with 
a convection-diffusion equation. However, perhaps because the strong consistency condition is 
routinely imposed in the construction of a numerical scheme, there are very few schemes, e.g., 
the L/D-F scheme, that are not strongly consistent with the PDE to be solved. 

At this juncture, note that the term "consistency" as defined on p.44 of [4] represents a 
concept that involves the numerical scheme, the PDE, and the rule of mesh refinement (i.e., how 
At and Ax are related as At, Ax 0). A two-level scheme is said to be consistent with a PDE if 
the truncation error -> 0 as At, Ax -» 0 under a given rule of mesh refinement (Note: For multi- 
level schemes, consistency means more than truncation error 0. See p.175 in [4], However, in 
this paper, the above definition of consistency for two-level schemes will be extended to a multi- 
level scheme like the L/D-F scheme. This extended definition should not be confused with the 
more rigorous definition used in an equivalence theorem given on p.172 of [4].) Generally, one 
can not say that a scheme is consistent with a PDE without specifying the particular rule of mesh 
refinement. A scheme may be consistent with a particular PDE under one rule of mesh 
refinement, and be consistent with another PDE under another rule of refinement. If a scheme is 
consistent with the same PDE regardless of how At, Ax -» 0, then it is strongly consistent with 
this particular PDE. 

As will be shown, the current numerical scheme is not a strongly consistent scheme. To expel 
any misconception that somehow such a scheme is intrinsically inferior than a strongly consistent 
scheme, next we shall compare the consistency, stability, convergence, and truncation errors of 
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two model schemes, i.e., the MacCormack scheme, which is strongly consistent with Eq. (2.2), 
and the L/D-F scheme, which is not. 

Let u(x,t) be a smooth function and uj u(JAx,nAt). Then with the aid of Taylor’s 
formula with remainder, it may be shown that 


and 


where 


[ FDE(M) ]" - [ PDE ]" = [ ER(M) ]" (ft > 0) 

[ FDE(L ) )] - [ PDE ] J = [ ER(L) ] J (\l 2: 0 ) 


( 6 . 1 ) 


( 6 . 2 ) 


[ PDE ] " & 


[FDE(M)]" W JL 


du du d 2 U 

n 

def 

du du d 2 U 

3 7 + a dx - 4 dx 2 

j 

dt dx ^ dx 2 


(6.3) 


J x = j Ax , / = nAl 

.«+ 1 5,5 s~n 1 , 5 2 5 t 5 2 v -» 

“/ - g ( 4 +,) “'- 2 + t + 


, . 5 35 2 2 \ 1/5 , 2 . 5x 5 2 

(1 - 1 + -r-" )u ' - 



5,5 

~ j( 4 -0«J« J 

(6.4) 

[ER(M)]" = <9 (At) + 0 [(Ax) 2 ] 

(6.5) 

[PDEtL)]- « ^ 

(1 + y )“/ - ( T + y)“/-i 



+ (x-|)i? + i - (l-y)«; '] 

(6.6) 


[ER(L)j; 


ft 


1 2 
At 

Ax 



_Af_ 

Ax 


2 

0[(Ar) 2 ] + 0[(Ar) 2 ] + <9 [(Ax) 2 ] 


V. J 


Several comments may be made about Eqs. (6.1) - (6.7): 


(6.7) 


a. 


[PDE]" = 0 (6.8) 

if u = u(x,t ) is a solution of Eq. (2.2). 
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b. u” = uj will satisfy the MacCormack scheme Eq. (D.3) if 

f FDE(M) ]" = 0 (6.9) 

As a result, [ FDE(M) ]" may be considered as the approximation of [ PDE ]" associated 
with the MacCormack scheme. Eq. (6.1) then states that [ ER(M) ]" is the error of this 
approximation. 

c. u] = m" will satisfy the L/D-F scheme Eq. (3. 19) if 

f FDE(L) ]" = 0 (6.10) 

As a result, [FDE(L )]" may be considered as the approximation of [PDE]" associated 
with the L/D-F scheme. Eq. (6.2) then states that [ ER(L) }* is the error of this 
approximation. 

d. Let u = u(x,t ) be a solution of Eq. (2.2). Then, by definition, [ FDE(M) ]] and [ FDE(L) ]] , 
respectively, are the truncation errors of the MacCormack and the L/D-F schemes at (j,n) 
[p.20, 4], Since Eq. (6.8) is satisfied, Eqs. (6.1) and (6.2) imply that 

[ FDE(M) ]" = [ ER(M) ]" (6.11) 

and 

[ FDE(L) ]" = [ER(L)j; (6.12) 

According to Eq. (6.5), [ ER(M) ]" — * 0 as At, Ax — > 0. Thus the truncation error 
[ FDE(M) ]" — > 0 as At, Ax — » 0. In other words, the MacCormack scheme is strongly 
consistent with Eq. (2.2). On the other hand, according to Eq. (6.7), [ER(L) ]", and thus 
the truncation error [ FDE(L) ] " , generally does not approach zero as At, Ax -» 0. As a 
result, the L/D-F scheme is not strongly consistent with Eq. (2.2). However, [ ER(L) ]J — » 
0 as At, Ax, At/Ax —> 0. Thus L/D-F scheme is consistent with Eq. (2.2) if the mesh is 
refined in a way such that At/ Ax 0 as Ar, Ax -> 0. 

e. Because (i) a discrete solution to Eq. (2.2) generated by an explicit scheme like the 
MacCormack scheme can not converge to its analytical counterpart without imposing the 
condition that At/Ax — > 0 as Ar, Ax — » 0, (ii) MacCormack scheme is strongly consistent 
with Eq. (2.2), and (iii) Lax’s equivalence theorem [p.45, 4] states that, given a properly 
posed initial-value problem and a finite-difference approximation to it that satisfies the 
consistency condition, stability is the necessary and sufficient condition for convergence, 
one comes to the conclusion that, in solving a properly posed convection-diffusion problem 
like that depicted in Fig. 6.2, stability of the MacCormack scheme will require that At/Ax 

0 os At, Ax — » 0. Note that the term "stability" referred to in Lax’s equivalence theorem 
has a meaning different from that defined in Section 4. 
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f. As will be shown in Appendix E, the Lax stability of the MacCormack scheme requires 
that the mesh be refined in a way such that the parameter 5 remains bounded as At, Ax -* 0. 
In the case where p > 0, this implies that 

■ = 0(Ax) (6.13) 

Ax 

It follows from Eq. (6. 1 3) that At /Ax — » 0 as At, Ax — > 0. Assuming that u = u(x,t ) satisfies 
Eq. (2.2), then Eqs. (6.5), (6.7), and (6.1 1) - (6.13) imply that 

[FDE(M)]" = <9 [(Ax) 2 ] and [ FDE(L) ]" = 0[(Ax) 2 ] (6.14) 

Thus the MacCormack scheme has no advantage over the L/D-F scheme in the order of the 
truncation error if the Lax stability is considered. 

This completes the comparisons between the MacCormack and the L/D-F schemes. It has 
been shown that a scheme that is strongly consistent may not have an intrinsic advantage over a 
scheme that is not. 


Next we shall study the consistency and the truncation error of the current scheme, 
particularly the two-level, two-dependent-variable discrete equations Eqs. (2.56) and (2.57). It 
will be shown that these two equations are consistent with a pair of PDEs if At l Ax — » 0 as At, Ax 
— » 0. One of these PDE’s is Eq. (2.2). 

To proceed, let u{x,t) and v(x,r) be smooth functions. Let 


- def - du 
W — V — — 

OX 

Letw" d - u(JAx,nAt), v" d — v(JAx,nAt) and w" d - w(JAx,nAt). Also let 


(6.15) 


. def _n+ 1 1 

[FI]] *2 Uj - ~ 


X + 


5(1 -T) 
1 - x 2 + 6 


,, s ~n l-X 2 -5 A ~n 

(1+X)U 7 _, + Ax Vj—l 


1 — X 1 


1 -x 2 + 5 


1 

+ 2 


X - 


(1 -x 2 )u, 

8(1 +x) 

1 -x 2 + 6 


x ( 1 - x 2 - 8 ) 


Ax v, 




( 1 — X) Uj+i - 


1 - x 2 - 8 


Atv, + , 


(6.16) 


» def _n+l 1 

[F 2]J =? A xv j + j- 


x + 


5(1 — T) 
1 -x 2 + 8 


1-x 2 


1 -x 2 -8 
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73e7s [ 4x “' + (1_x2 “ 5)Ax ^] 


1 - x 2 - 5 


1 5(l+x) „ 1 — x 2 -« /1 , „ 1-x 2 -5 a 

+ l[ T_ 7^771 ' <1+X) 1T775 J>l 

Then, with the aid of Taylor’s formula with remainder and the assumptions that 


It may be shown that 


5 > 0 and 1 > x 2 


[FDE1 ]" - [ PDE ]" s [ER1 ]" 


[ FDE2 ]J - w" = [ER2 ]" 
Here [PDE ]" is defined in Eq. (6.3). Also 

[FDE1 ]" & ~ [F 1 ]] 


[ FDE2 ]" % - (1 t2 + 5)2 [F2]] 
1 45(l-x 2 )Ax j 


[erhj | * 

1 K 1 - x 2 + 8 dx 2 dt 2 dx 2 

V. J j ^ ' J ^ J J 


4 dx 2 


+ f \m\ (Ax) 2 -11^1 -C1 ^ t 2 -5)[(Ax) 2 -^ 2 ( Ar) 2 ] [ + 0[m2 . 

O dx J 
^ J J 


1 -x 2 + 5 


+ O [ (Af )^ ] 


+ X + + g -^-|(l+x)xO[(Ax) 4 ] + (1 -x 2 -5)x0[(Ax) 4 ]J 


+ T ~ f-T^+S (l-x)x0[(Ax) 4 ] + (1 -x 2 -5)x0[(Ax) 4 ] I (6.23) 
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def (1 -T 2 + 5) 2 (Ax) 2 

r ”> 

o(l - x 2 -5) 2 (Ax) 2 

> 

OJ 

/ 

16|i(l -T 2 ) 

dt 

» j 

16p(l -T 2 ) 
j 

l dx ) 


, a( 1 ~t 2 )(Ax) 2 
8|x 



dx 2 . 


+ 


( 1 - T 2 + 8 ) 2 

4 8 ( 1 - T 2 ) 


x O [ (Af) 2 ] 


+ 


+ 


1 + f ci-" 2 ) 

l-f(l-T’) 


x jo [(Ax) 2 ] + 

1 -T 2 -8 
1 + T 

x 0 [ (Ax) 2 ] | 

x \ 0 [( Ax) 2 ] + 

1 - T 2 - 8 ' 

1 -x 

k. J 

x 0 [ (Ax) 2 ] • 

J l 

j 


(6.24) 


The significance of Eqs. (6. 19) - (6.24) will be discussed in the following comments: 


a. [FDE1 ]" and [FDE2 ]", respectively, may be considered approximations of [PDE]" and 
w". Eqs. (6.19) and (6.20) then state that [ER1 ]" and [ER2]", respectively, are the errors 
of these approximations. 


b. Eqs. (2.56) and (2.57), respectively, are equivalent to 

[FDEl ]" = 0 (6.25) 

and 


[ FDE2 ]" = 0 (6.26) 

if u" and v" in Eqs. (6.25) and (6.26), respectively, are replaced by and a". 

c. With the aid of the observations made in (a) and (b), and Eqs. (6.3) and (6.15), one 
concludes that Eqs. (2.56) and (2.57), respectively, may be considered as the discrete 
approximations of Eq. (2.2) and the PDE 


v 



(6.27) 


with the understanding that and a", respectively, are the discrete counterparts of u and v. 

d. Let u = u(x,t) and v = v(x,r) be a solution of the system of PDEs Eqs. (2.2) and (6.27). 
Then, by definition, [FDEl]" and [FDE2]" are the truncation errors of the discrete 
equations Eqs. (2.56) and (2.57). Furthermore, we have [ PDE ]" = 0 and w" = 0. Thus 

[FDEl ]" = [ER1 ]" and [ FDE2 ]" = [ER2]" (6.28) 
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i.e., [ER1 ]" and [ER2 ]" become the truncation errors. 

e. Since p and a are constants, parameters t and 5 vary as At and Ax vary. Generally, z and 8 
do not approach certain limits as At, Ax — » 0. In spite of this and the fact that [ER1 ]" and 
[ER2 ]" are dependent on x and 5, it is shown in Appendix E that 

n.St,o [ER11 ’ = 0 <«•*» 

and 


lim [ER2 ]7 = 0 

A//Az , Ax — > 0 1 


(6.30) 


if u = u(x,t) and v = v(x,r) satisfy Eq. (6.27). Note that the notation At -» 0 does not 
appear on the left side of Eq. (6.30) because At — > 0 if At/ Ax, Ax — > 0. Furthermore, it is 
shown in Appendix E that 


[ERl]; = 0[(Ax) 2 ] and [ ER2 ]* = O [ (Ax) 2 ] (6.31) 

if u = u(x,t ) and v = v(x,r) satisfy Eq. (6.27) and the rule of refinement is such that 6 
remains bounded as At, Ax — > 0. 


f. Combining (a) - (e), one may conclude that Eqs. (2.56) and (2.57) are consistent with Eqs. 
(2.2) and (6.27) if the rule of mesh refinement is such that At/Ax — » 0 as Ar, Ax — » 0. 

g. Because of Eq. (6.27), convergence of a discrete solution to its analytical 

counterpart ( u , v ) implies that y J u and ot™ — > du/dx. This is consistent with the 

interpretation ofy? and a] given in Eqs. (1.17a) and (1.17b). 

Finally we consider the special case in which 8 = 0 (i.e., p = 0). For this case, it may be 
shown that 


A7 tF11 '" - 


3m , 3m 

n 

At 

3 2 m 

2 3 2 m 

— + a — 

•>' Sx i 

j 

2 

3r 2 

3x 2 

•4 


+ 0[(Ax) 2 ] + 0[(At) 2 ] (6.32) 


1 r/79 _ 

dw 

3vv 

. A 

' 3m 

, n 3m 

Af At [F1] ' 

3r 

a 3x 

dx 

dt 

3x 

J . 


A i 
2 


3 2 vv 

2 3 2 w 

3 

r d 2 u 

2 3 2 m 

dt 2 

3x 2 

dx 

[ 3r 2 

3x 2 J 


+ O [(Ax) 2 ] + O [ (Ar) 2 ] (6.33) 


Using Eqs. (6.32) and (6.33) and some comments given previously, one concludes that Eqs. 
(2.56) and (2.57) are strongly consistent with Eq. (1.3) and 



du du 

— + a^— 

dt dx 


(6.34) 


3w 3w d_ 

dt 0 dx dx 



Eqs. (1.3) and (6.34) are equivalent to Eq. (1.3) and 


dw 

dt 



(6.35) 


Eqs. (1.3) and (6.35) are similar in their forms. They, respectively, represent the wave motions 
propagating with the speeds +a and -a. 

Let u = u(x,t) and v = v(x ,t) satisfy Eqs. (1.3) and (6.35). Then the lowest-order terms on the 
right sides of Eqs. (6.32) and (6.33) vanish. Thus Eqs. (2.56) and (2.57) are strongly consistent 
with Eqs. (1.3) and (6.35) with their truncation errors being O [( Ax ) 2 ] + O [ (At ) 2 ]. 


- 63 - 



7. NUMERICAL EVALUATION 


In this section, the current method will be compared numerically with the L/D-F method and 
the MacCormack method. During this comparison, many theoretical results developed 
previously will also be evaluated using the numerical results presented. 

To simplify our effort, the numerical problems to be considered have the common initial 
condition 

u° = sin(2njAx) , j = 0, ±1, ±2, • • • (7.1) 

where 



(7.2) 


and K > 3 is the integer introduced in Eq. (4.1). According to Eqs. (7.1) and (7.2), u° is periodic 
over every unit length and every K mesh intervals in the Jt-direction. In this section, we shall 
continue to use the mesh depicted in Fig. 2.1(a) and assume the periodic condition, i.e., u" = 
u] +K ,j = 0, ±1, ±2, ••■,/! = 0, 1,2, •••. 

The analytical problem corresponding to the above numerical problem will be specified by 
the periodic condition u(x+KAx,t ) = u(x,t) and the initial condition u(x, 0) = I(x) where f(x) is 
determined according to Eqs. (5.7) and (5.8). With the aid of Eqs. (4.7), (5.7), (5.8), (7.1), and 
(7.2), it may be shown that 


t>Q — 0 , b ] 



b 2 - bj ~ ■ = bf(-2 = 0 » 


b K -\ = ~ 


2 i 



0a:-i = - 


2jc 

K 


(7.3) 

(7.4) 


and 


I(x) = sin(27tx) (7.5) 

By using Eqs. (5.15) and (7.2) - (7.4), one obtains the analytical solution to Eq. (2.2), i.e., 

u = u a (x,t ) =? sin[2jc(^-a/)] (7.6) 

Note that parameters yj and a°, which are needed in the initiation of the current numerical 
procedure, may also be evaluated by using Eqs. (5.1), (5.4), (7.1), and (7.5). 

Combining Eqs. (5.17), (5.52), (5.53), and (7.2) - (7.4), it may be shown that 

u a {x n j,t n ) = e" 4 ^ nA/ sin{2rt[yAx-n(a-fe)Ar]) = u a (j,n, 1) + u a (j,n,K-l) (7.7) 

with 
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I 11 


(7.8) 


u a (j,n, 1) = jjD(nAt)e iri 


and 


u a (j,n,K-\) = --^D(nAt)e ^ (7.9) 

2i 

Here 

D(t ) e~ 4 ^' (t> 0) (7.10) 

is the decay factor of u a (x,t ) and 

V ^ 2%{j-nx)IK (7.11) 

is a phase angle. u a (xj,t n ) is the value of the analytical solution at the mesh point (J,n). Its 
numerical counterpart in the present method is f-. By using Eqs. (2.23), (2.24), (5.13), (5.14), 
(5.54), and (7.3), one concludes that 


Yy = «(/>, \) + u(j,n,K-\) 


(7.12) 


Let 


„ , def ij-U a (x n j,t n ) 

R\ (/.«) - ~ tt — 

D ( nAt ) 


(7.13) 


In other words, R \(j,n) is the error of the numerical solution f] normalized by the decay factor 
ofu a (x n j,t n ). Eqs. (5.55) and (7.7) - (7.13) imply that 


*.(/.«) = 77 

2 / 


r x {n,\)e lt>1 -r x {n,K-\)e 


(7.14) 


Substituting Eqs. (5.46), (5.48), and (5.49) into Eq. (7.14) and using Eq. (7.4) one has 
' T Ai(x,8) 


Riij.n) = n- 


192 


24 

1-5t 2 

l-T 2 


(2tc/AT) j cos (<!> ; n ) 


Ai(t,5)-4t 2 


( 2jc / K ) 4 sin ( <t>" ) + O [ (2jc /K ) 5 ] 


(1-T 2 ) 2 -S 2 j x8 „ ivs 3 

i x-(2k K) cos(<b") 

32( 1 -T 2 ) 1-T 2 


_ 1 _ 

24 


( 1 + 3 T Z ) + 


2 \ . 3 5 2 (1-9 t 2 ) 


d-T Z ) 


2 \2 


(2% IK) 4 sm(^1) + 0[(2n/K) 5 ] 


(7.15) 
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According to Eq. (5.56), u ax (x],t n ) is the value of du a (x,t)/dx at the mesh point (j,n). Its 
numerical counterpart in the current method is a". By using Eqs. (5.17), (5.56), (5.57), and (7.2) 
- (7.4), one concludes that 

Uax(x],t n ) = 2ne- A * 2 » n * cos[2n[jAx-n(a-b)&t]) = u ax (j,n,\) + u^lj^K-X) (7.16) 
with 


»ax(j\n, 1) = nD(nAt)e‘^ 
Uax(j,n,K-\) = nD(nAt)e~‘* j 


Moreover, Eqs. (2.23), (2.24), (5.13), (5.14), (5.58), and (7.3) imply that 


CXy — {£*(/, w, 1) + u x (j,n,K 1) 


Let 


(7.17) 

(7.18) 


(7.19) 




def 


a" ~ UgxOCjJ^) 
2nD(n At) 


(7.20) 


i.e., R 2 (j,n) is the error of a” normalized by the "amplitude" of u ax (x” ,t n ). Then Eqs. (5.59) and 
(7.16) - (7.20) may be used to show that 


*2 (/.«) = J 


r 2 (n, l)e l<t>; +r 2 (n,K- l)e 


-‘*7 ' 


(7.21) 


Substituting Eqs. (5.46), (5.50), and (5.51) into Eq. (7.21) and using Eq. (7.4), one has 
tAi(t,5) 


A 2 (j.n) = n|- 

8 

+ 

192 


24 

1-9t 2 
1 -t 2 


(2n/ K ) 3 sin ( <J)y ) 


A 1 (x,8)-4t : 


( 2k/K f cos ( + O [ ( 2n/K ) 5 ] 




48 


(1+3t 2 ) + - §2(1 j 
(1-X 2 ) 2 


( 2n /K ) 2 cos (<t>; ) + O [ ( 2n /K ) 3 ] 


(7.22) 


According to Eqs. (7.3), (B.31), and (B.33), the numerical counterpart of u a (x n j,t n ) in the 
L/D-F method is 




(7.23) 


Let 
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n . def !L(j,n)-u a (x],t n ) 

R L (j>n) - pT7 

D (nAt) 


Then Eqs. (5.60), (7.7) - (7.9), and (7.23) imply that 

1 


fl t (/.«) = 


2i 


(n, l)e*' -r L (n,K-\)e~^ 


(7.24) 


(7.25) 


Substituting Eqs. (5.61) - (5.63) into Eq. (7.25) and using Eq. (7.4), one has 
R L (j,n) = (2n I K ) 2 sin(<{>" ) + -^- ( 1 -x 2 )( 1 - -j 5 2 )(2tc/AT) 3 cos ( <J>" ) 


96 


95 2 x 2 (l-x 2 ) + 35x 4 +4x 2 -y 5 2 + 2 


( 2k! K ) 4 sin (<t>?) + O [( 2k/K ) 5 ] 


+ {d )j x 2 (2 k IK) 2 sin(<j>” ) - ( 1 — 2 x 2 ) ( 2 jc / A! - ) 3 cos ( <()" ) 


48 


4x 2 (4-9x 2 ) + 35 2 (15x 4 -12x 2 + l) j (27C//S: ) 4 sin(<(>" ) 


-t- O [ (2 tc/A' ) 5 


(7.26) 


According to Eqs. (7.3), (D.5), and (D.6), the numerical counterpart of u a (xj,t n ) in the 
MacCormack scheme is 


«*(/» = HM(j>n,\) + H M (j,n,K-\) 


Let 


(7.27) 


_ s def u M (J,n)-u a {x n j ,t n ) 

R*(j.n) ~ n ( ~\T\ 

D (nAt) 


Then Eqs. (5.64), (7.7) - (7.9), and (7.28) imply that 


r„(n, 


*•«•*>- 2 ?L 


Substimting Eq. (5.65) into Eq. (7.29) and using Eq. (7.4), one has 
R*(j \n) = /1 | ~ 1 6 T ~ ( 2n / ^) 3 cos (<(>") 


(7.28) 


(7.29) 
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+ — [ 5( 1 -6x 2 ) - 6x 2 ( 1 -x 2 ) ] (2n/K ) 4 sin(<t>; ) + O [(2jc/A ) 5 ] 


Let x = 0. Then Eqs. (7.15), (7.22), (7.26), and (7.30), respectively, are reduced to 


R\(j,n) = n\ 5(1 19 2 52) (27c/*) 4 sin(<l)") + 0[(2re//O 5 ] 


(l-5 2 )(l+35 2 ) 


(2re//sr ) 4 sin (*)>;) + (9 [(2tc/ /ST ) 3 


(x = 0) 


Ri(j,n) = J 5(1 ^ 52 ) (27t//:) 4 cos(4>?) + 6>[(2n/^) 5 


+ ( 2k 7 K > 2 cos ( ♦" ) + O [ ( I K ) 3 1 


(x = 0) 


R L (j,n) = «|^^|^-(27c/A) 4 sin((t);) + <9[(2jc/^) 5 ] 


^ ( 1 - Y )(2nfK ) 4 sin (*?) + <9 [ (2k/ K) 5 


(x = 0) 


R„(j.n) = n< -^(2K/K) 4 sm(^) + 0[(2n/Kf] 


(x = 0) (7.34) 


Let x and 5 satisfy the optimal condition Eq. (5.70), then Eqs. (7.15) and (7.22), respectively, 
are reduced to 


R\(j,n) = nj - T ~~j= ~ (2^/A) 4 sin(<j>") + 0[(27i/A') 5 


^ 4 %-\ 7^^ ( 2tc/ AT ) 3 cos(<|>" ) + -jy ( 1 - 3x 2 )(2jt/A r ) 4 sin ( <)>" ) 


+ <9[(2tc/A) 5 


(7.35) 
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R 2 (j,n) = n |- — ^ ( 2n I K ) 4 cos ( tf) + O [ ( 2k /K ) 5 ] 

- ^=r ( 2n/K ) sin (*J ) + -^- ( 1 - 1 2 ) ( 2k/K ) 2 cos (*? ) 

+ O [(27C/AT) 3 ] (l-t 2 = V3 5) (7.36) 

According to Eqs. (5.46), (5.48), and (5.49), r x (n,k) may be approximated by the sum of the 
expressions on the right sides of Eqs. (5.48) and (5.49). This sum may be converted to the 
approximate form of R \ ( j,n ) given in Eq. (7.15) if one carries out the following substitution: 

(2k/ K y sin ( <}>" ) if / is even 

(e*y -» \ 

— ( 2tc / AT y cos ( <t>7 ) if / is odd 
^ i 1 

The same relation also exists between r L (n,k ) and R L (J,n), and between r M (n,k) and R M (j,n). A 
similar relation also exists between r 2 (n,k ) and R 2 (j,n). However, for this case, the rule of 
substitution is 

f (2k /K) 1 cos(<j>" ) if l is even 1 

(0*)' -> i \ (7.38) 

[ i(2n/K) sin(0" ) if / is odd J 

The approximations for the spurious parts r \-(n,k), r 2 ~(n,k), and r L _(n,k), respectively, given 
by Eqs. (5.49), (5.51), and (5.63) are used in the above derivation of Eqs. (7.15), (7.22), and 
(7.26). In the case where 5 = 0 or 8 is very close to zero, these spurious-part approximations may 
no longer be valid. Fortunately, the spurious part generally is negligible compared with the 
principal part in a calculation with large n. Thus we may completely omit the spurious parts in 
Eqs. (7.15), (7.22), and (7.26) if 5 is very close to zero and n is large. 

Subject to the modifications required by the substitution rule given in Eqs. (7.37) and (7.38), 
the comments made in Section 5 about the approximate forms of r x (n,k), r 2 (n,k), r L (n,k), and 
rti(n,k) are applicable to those of R\(j, n), R 2 (j,n), R L (j,n), and R u (J,n). Particularly, the leading 
terms in the principal parts of R x ( j,n ) and R 2 (J,n) also vanish if A|(t,5) = 0. Because the factor 
2k/ K — > 0 as K — » +°° (i.e., Ax -4 0), the leading term becomes more dominant as K — > +«. As 
a result, the effect of the leading-term annihilation on the accuracy of a numerical method 
becomes more pronounced as AT —>+<». 

Eq. (7.6) implies that 
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du a (x,t) 

dx 


= 2 nae 4n ^' cos[2jc(j:-af)] 


and 


d 2 u a (x,t) _ . 2 

\i — — 2 — = -4x 2 n«- 4 *»“ sin[27i(x-a/)] 

OX 


(7.39) 


(7.40) 


Thus, in the case where u = u a (x,t), the relative importance of the convection and the diffusion 
terms in Eq. (2.2) may be determined by comparing the "amplitude" of the expressions on the 
right sides of Eqs. (7.39) and (7.40). As a result, in the following numerical study, the solution 
u = u a (x, f) will be referred to as: (i) convection dominant if |<z| 2 7 i p, ( ii) convection- 

diffusion comparable if \a \ ~ 2 rcp, and (iii) diffusion dominant if \a | c 27 tp. 

In the following numerical study, each test problem is defined by the values of (i) the physical 

parameters a and p, and (ii) the mesh parameters b, n. At, and K (= — ). After n time steps the 

Ax 

numerical error of a test problem will be measured by 


E(b,n,At,K) & logic 


AT-I 


[ D(nAt) K 


Y Z \u]-u a (x],t n )\ 


(7.41) 


where u" is the numerical solution at the mesh point (j,n). Roughly speaking, the negative of 
E(b,n,At,K ) represents the average number of correct significant figures in u", 
j — 0, I, 2, • ■ • , K— 1. Note that, in the current paper, we will not distinguish between the exact 
solution of a numerical scheme and the actual numerical solution of the same scheme, i.e., the 
roundoff error is assumed to be negligible. 

With the aid of Eqs. (7.13), (7.24), and (7.28), one concludes that 

log to R l ( b , n, A t,K) ( Present scheme ) 


E(b,n,At,K) = \ logio R L (b,n,At,K) ( L/D-F scheme ) 


logio R\t(b,n, A/,*) ( MacCormack scheme ) 


(7.42) 


Here 


R \(b,n,At,K) V I«t0‘.«)| 

K j-Q 

R L (b,n,At,K) i ¥ -i- V |/?l0‘.«)| 
K j = o 


(7.43) 

(7.44) 


and 
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(7.45) 


J 

R u {b,n,At,K) & ±r £ !**(/.«) I 

A j=0 

, respectively, are the averages of | /?,(/, /i)|, \R L (J> n )\ and | /?*(/, n)| at time level n. Note that 
R i (j,n), R L (j,n), and R u (j,n), implicitly, are functions of b. At, and K. In our future discussions, 
E{b,n,At,K), R x (b,n,At,K), R L (jb,n,At,K), and R M (b,n,At,K), respectively, may be abbreviated 
as E, R i , R l , and R„. 


To pave the way for the interpretation of the numerical results to be presented, several 
approximations of R i , R L , and R u will be introduced in the following discussions: 

a. Let 1 - t 2 = VT 5. Then the approximation of Ri(j,n) given in Eq. (7.35) may be 
applicable. Let R \ ( j,n ) be dominant by the leading term in its principal part. Then it may 
be shown that 


- ^ 2 n jc 3 x 2 1 1 - x 2 1 


(1 -t 2 =V3 5 , n, K» 1) 


(7.46) 


Note that, in obtaining Eq. (7.46), we use the identity 


1 


K - 1 


„ lim \ Y s |sin(<t>; + C)| 

k -»+oo y k j=0 


1_ 

K 


(7.47) 


where C is any real number which is independent of j. 


The parameters Re, To, 5 q, and (Ar) 0 were defined in terms of (a -b), p, and Ax 
(= 1 IK) in Section 5. By these definitions, we have 1 -Xq=^ S 0 . Furthermore, Eqs. 
(5.73) - (5.75) imply that 


|x 0 | < 1 and 5 0 = -4=- if |Re| <e 1 


(7.48) 


Combining Eqs. (2.17), (2.18), (7.46), and (7.48), one has 

( n, K » 1 » | Re | ) 


d th r\t\ • ^(a-bfnmo 
Ri(b,n,(At) 0 ,K) = - 


18 p IT 


(7.49) 


i.e., R\(b,n, {At) 0 ,K) is approximately proportional to n (Ar)o (i.e., the total running time), 
and (Ax) 4 if n, K >• 1 » j Re j . 

b. Let t 2 < 1. Then with the aid of Eqs. (2.17), (2.18), and (7.2), Eq. (7.30) may be further 
simplified as 


= Va 2 + A 2 sin ( +£q) + n 0[(2n/K) 5 ] (t 2 «1) (7.50) 


where 
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(7.51) 


4 <kf 47C 3 (a-6)(« A/) 4 def 4 Jt 4 p. ( n Af ) 

A 1 - 7773 * a 2 7775 


3 K 


3K' 


and £o is a number such that 


A i A 2 

smeo = - ,_ A . „ , - , cos eo = 


Va? +a| 


Va? +a! 


(7.52) 


Assuming that the term n O [(2n/K ) 5 ] in Eq. (7.50) is negligible, then Eqs. (7.45), (7.47), 
(7.50), and (7.51) imply that 


» A 8re 2 V( a -^)2 + ^ 4 2 (nAr) 

- 7^2 — (AT»l»r) 


(7.53) 


i.e., R m is proportional to nAt and (Ax) 2 if K » 1 » x 2 . 


c. Let R\{j,n ) given in Eq. (7.15) be dominated by the leading term in its principal part. 
Then it may be shown that 


R i 


2n n 2 |tAj(t,5)| 

3K 3 


(/ST»1) 


(7.54) 


Similarly, let R u (J,n) given in Eq. (7.30) be dominated by its leading term. Then it may be 
shown that 


Rm = 


8»7t 2 I t ( 1 — t 2 ) I 
3 a : 3 


(tf 1) 


(7.55) 


Let n x 0. Then the right sides of Eqs. (7.54) and (7.55) will be equal if and only if either 
(i) ( 1 - x 2 ) 2 + 5 2 = 0 or (ii) 5 ( 1 - x 2 ) 2 - 3 5 2 = 0. Since 1 - x 2 > 0 and 8 > 0, case (i) 
cannot occur and case (ii) is reduced to 

1 - x 2 = ^3/5 8 (7.56) 

i.e., R 1 and R„ are approximately equal if x and 5 are related by Eq. (7.56). Note that the 
only difference between Eqs. (5.70) and (7.56) is that the factor ^3 in the former is 
replaced by ^3/5 in the latter. 


Let ( a - b ), p, and Ax be given. Then it may be shown that there is one and only one 
set of values of x, 5, and At (denoted by Xq, 8q, and (A/)o, respectively) such that Eq. (7.56) 
is satisfied. These values may be determined by using a procedure similar to that used in 
the determination of x 0 , 80, and (Ar)o (see Section 5). 

Let p, Ar, and Ax be given. Eq. (7.56) has no solution if 8 > V5/3. if V5/3 ^ 5 > 0, 
then Eq. (7.56) is satisfied if 


TT 
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X = X± *=f ±Vl _ V 3/5 5 


(7.57) 


d. Let x = 0. Then the approximations given in Eqs. (7.31), (7.33), and (7.34) may be 
applicable. In each of them, we further assume that it is dominated by the leading term in 
its principal part. Then it may be shown that 


- . h ft 3 5 j 1 ~ 3 1 , 2^n(„Arj | l-36 z | 

6 K 4 3 K 2 


Rl = 


2 n k 3 5 1 1 - 8 2 1 8 7t 3 p ( n At ) 1 1 - ^ 8 2 | 


3 A" 4 3 AT 2 

- ^ 2nn 3 5 8rc 3 p(rtAr) 

M ” 3 AT 4 " 3 K 2 


(x = 0, AT » 1) 


(x = 0, A' 3> 1 ) 


(7.58) 


(7.59) 

(7.60) 


Thus R m , approximately, is proportional to nAt and (Ax) 2 if x = 0. However, one must 
assume both x = 0 and 1 5| 1 in order that R 1 and R L , be approximately proportional to 

nAt and (Ax) 2 . 


Assuming 5 > 0, one may conclude that 


(0 

>11 

1* 

if 5 = V 5 / 6 =0.913 

(x = 0, AT» 1 ) 

(7.61) 

(«) 

Rl = Rm 

if 5 = ^5/J = 1.291 

(x = 0, AT » 1) 

(7.62) 

and 





m 

r*' 

II- 

if 5 = ^873 = 1.633 

(x = 0, AT » 1 ) 

(7.63) 


This completes the preliminaries for the following numerical evaluation. Several sets of test 
problems will be considered. In set #1, the test problems have the same values of a (=1), p 
(=0.1), b (=0), K (=30), and t (=0.5). They are distinguished by their values of At. For each 
member of set #1 , n = IN(r /At) where 

def 

IN(x) - the integer nearest to the real number x. (7.64) 

As a result, nAt = t if t /At » 1. Sets #2 - #4 are defined similarly. The values of a, p, b, K, t, 
Re, Xo, x,, and Xq for sets #1 - #4 are listed in Table 7.1. Among these parameters, x* is the only 
one yet to be defined. 

For sets #1 - #3, | | = 2jxp. Thus the corresponding u a (x,t ) is convection-diffusion 
comparable. For set #4, | a \ » 2rtp. Thus the corresponding u a (x,t ) is convection dominant. 

def 

Because (i) x - (a - b ) At /Ax, and (ii) a, b, and Ax are constant among the test problems in 
any one of sets #1 - #4, each test problem within one set can be identified by its value of X. In 
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Figs. 7.1 - 7.4, the values of E(b, IN(r /At),At,K) for the current, the L/D-F, and the MacCormack 
schemes are plotted against the values of x for representative test problems in sets #1 - #4. The 
following remarks pertain to the results shown in these figures: 

a. For the test problems in sets #1 and #2, 



Thus, on the 5-x plane, each of these problems is represented by a point on the straight 
line x = 5/12. Since [Re | 1, one may conclude from Fig. 4.1 that the MacCormack 

scheme will become unstable at a x with | x | «c 1. Asa matter of fact, it may be estimated 
that the above straight line intercepts the stability boundary at x = 0.17. Note that, as a 
result of how stability is defined (see Section 4), instability generally occurs at a value of x 
slightly greater than what is predicted by the stability map \iK and n are finite. 

For the test problems in set #3, Re = 1/24. Thus the MacCormack scheme becomes 
unstable at an even smaller x. On the other hand, Re = 5/6 = 1 for the test problems in set 
#4. Thus the MacCormack scheme becomes unstable at a x very close to I. 

b. The most remarkable result shown in each of Figs. 7.1 - 7.4 is the sharp increase in the 

accuracy of the current method in the neighborhood of x = x 0 . For each of sets #1 - #3, 
which is characterized by |Re | 1, the peak accuracy is 3 - 4 orders of magnitude higher 

than what can be achieved by either the L/D-F scheme or the MacCormack scheme. For 
set #4, which is characterized by Re = 1, the peak accuracy of the current method is about a 
factor of 30 higher than that of the MacCormack method which occurs just before it 
becomes unstable. 

Moreover, for each of sets #1 - #4, the actual value of x (denoted by x* in Table 7.1) at 
which the peak accuracy of the current method occurs is extremely close to the theoretical 
value x 0 . The difference is numerically insignificant. As an example, it can be determined 
numerically that Eq = —4.923 and E, = —4.924 for set #1. Here Eq and E t , respectively, 
denote the values of E of the present method for the test problems with x = x 0 and x = x, . 

c. For sets #1 - #3, | Re | c 1. Thus Eq may be estimated by using Eqs. (7.42) and (7.49). 
Using these equations, one obtains E 0 = -4.97 for set #1, E 0 = -4.67 for set #2, and E 0 = 
-6.18 for set #3. For set #4, | Re | = 1 and thus Eq. (7.49) is not applicable. However, by 
using Eqs. (7.42) and (7.46), one has Eq = —3.33 for set #4. From Figs. 7.1 - 7.4, one can 
infer that the above estimates of E 0 agree very well with those of £*. 

d. According to Eqs. (7.42) and (7.53), parameter E obtained using the MacCormack scheme 
is not dependent on At (and thus x) if ( a - b ), p, K, and nAt are held constant and K » 1 ;*> 


HI 
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x 2 . This explains the fact that, for the MacCormack scheme, the values of E shown in each 
of Figs. 7. 1 - 7.3 are hardly dependent on x. 

Moreover, the approximation of R M given by Eq. (7.53) generally is very accurate. As 
an example, for the test problems in sets #1 - #3 with x = 0.01, estimates of E obtained by 
using Eqs. (7.42) and (7.53) are -1.8146, —1.5135, and —2.4166, respectively. On the other 
hand, the actual numerical values are -1.8148, -1.5119, and -2.4167, respectively. The 
above comparison also indicates that the accuracy of the approximation of R M improves 
with the absolute value of E, i.e., the accuracy of the numerical calculation. The reason for 
this fact is explained in statements following Eqs. (5.40) and (D.7). 

e. It was argued in Section 5 that the current scheme will generally be more accurate than the 
MacCormack scheme by a factor of 4 when 5 is small and the initial condition is smooth. 
Since x = Re-5 and Re is a constant in each of sets #1 - #4, one would expect that, as x — > 0, 
the value of |£| for the current scheme will be greater than that for the MacCormack 
scheme by approximately logio4 = 0.602. This expectation is certainly confirmed by the 
results shown in Figs. 7.1 - 7.4. The same results also indicate that, as x — > 0, the value of 
E for the MacCormack scheme coincides with that for the L/D-F scheme. This fact may be 
explained by the following observations: The first term in the principal part of R L (j,n) 
given in Eq. (7.26) becomes negligible compared with the second term as x — » 0 (and 5 = 
x/Re 0). On the other hand, the latter is reduced to the leading term in R u (j,n ) (see Eq. 
(7.30)) as 5 — » 0. 

f. The values of Xo for sets #1 - #4 are given in Table 7.1. From Figs. 7.1 - 7.4, it is seen 
that, for each of sets #1 - #4, the value of E for the current scheme is approximately equal 
to that for the MacCormack scheme at x = Xo. This observation confirms a prediction made 
in a previous theoretical discussion. 

This completes the discussion of the test problems in sets #1 - #4. Next we will study the test 
problems in sets #5 and #6. As shown in Table 7.2, the test problems in each of these sets share 
the same values of a, p, K, n, and At. They are different in their values of b. As a result, each 
member in set #5 or #6 may be identified by its value of x. Note that, for each test problem in 
sets #5 and #6, the total running time nAt = 26/(3^) = 5.0037. 

Since a-b = 1.0 for sets #1 - #4, no member in these sets may have x = 0. Contrarily, 
(a - b ) is a variable in sets #5 and #6 and therefore both contain a member with x =0. Recall 
that the leading terms in the principal and spurious parts of R x (j,n), R 2 (J,n),R L (J,n) and R M (J,n) 
will be annihilated if x = 0. In other words, for either of sets #5 and #6, and for all three schemes 
considered, the test problem with x = 0 is expected to have the highest accuracy. 


- 75 - 



The leading terms in the principal parts of Ri(j,n) and R 2 (J,n ) can also be annihilated if t = 
T + or x = x_. The parameters x + and x_ are defined in Eq. (5.77) and the values of x + for sets #5 
and #6 are given in Table 7.2. Note that, for set #6, x + = x_ = 0. Thus the accuracy peaks of the 
current scheme that occur at x = 0, x = x + , and x = x_ will merge into one. Furthermore, when x = 
x + = x_ = 0, Eq. (7.15) is reduced to 

fliO'.rt) = n O [(2n/K) 5 ] - ~—(2n/K) 4 sin(<j>J ) + 0 [(27 c//sT ) 5 ] (7.65) 

i.e., all the explicitly-given leading terms in the principal part vanish. As a result, the 
approximation of R i given by Eq. (7.58) also vanishes. Obviously, in the case where x+ = x_ = 0, 
the estimation of E at x = 0 for the present scheme requires an approximation containing more 
explicitly-given terms. 

In Figs. 7.5 and 7.6, the values of E for the current, the L/D-F, and the MacCoimack schemes 
are plotted against the values of x for representative test problems in set #5 and #6. The 
following remarks are for the results shown in these figures: 

a. On the 5-x plane, each test problem in set #5 is represented by a point on the vertical line 
5 = O.S/^3. From Fig. 4.1, one concludes that the MacCoimack scheme will be stable for 
a test problem as long as |x| < 0.9. A similar conclusion may be applied to the test 
problems in set #6. These observations are consistent with the numerical results shown in 
Figs. 7.5 and 7.6. 

b. As expected, the accuracy of all three schemes considered reaches its highest level at x = 0. 

The values of E at x = 0 can be estimated by using Eqs. (7.42) and (7.58) - (7.60). For 
the test problem with x = 0 in set #5, we obtain E = -3.38, -2.41, and -2.34, respectively, 
for the current, the L/D-F, and the MacCoimack schemes. These estimates agree very well 
with the results shown in Fig. 7.5. 

For the test problem with x = 0 in set #6, Eqs. (7.42), (7.59), and (7.60) imply that E = 
-2.46 and -2.34, respectively, for the L/D-F and the MacCoimack schemes. These 
estimates are also in good agreement with the results shown in Fig. 7.6. As noted 
previously, for the test problem considered here, the approximation given in this paper is 
inadequate in providing an estimate of E for the current scheme. 

c. In Fig. 7.5, another accuracy peak for the current scheme also appears in the neighborhood 
of x = x + . The actual value of x (denoted by x, in Table 7.2) at which the peak accuracy 
occurs again is very close to the theoretical value x + . With the aid of Eqs. (5.77), (7.46), 
and (7.42), it is estimated that E = -3.04 for the present scheme at x = x+. This is very 
close to the actual peak value one observes in Fig. 7.5. Note that the accuracy peaks at x = 
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x + and x = 0 are merged in Fig. 7.6. 

d. From Figs. 7.5 and 7.6, and Table 7.2, one may confirm that, for each of sets #5 and #6, the 
values of E for the current and the MacCormack schemes are approximately equal at x = 

V 

Next we study the test problems in set #7. According to Table 7.3, these problems share the 
same values of a , |i, b, K, and t. They differ in their values of At. Again we assume that n = 
IN(r/Ar). Because the relations x = 0 and 5 = 36 Ar hold for these problems, the error measure E 
is plotted against 5 in Fig. 7.7 for all three schemes considered. From this figure, one observes 
that 

a. The value of E for the MacCormack scheme is hardly dependent on 8 before it becomes 
unstable near 8 = 2.1. 

b. As 8 —> 0, (i) the value of E for the L/D-F scheme approaches that for the MacCormack 
scheme, and (ii) the difference between the value of E for the current scheme and that for 
the other two schemes approaches logi 0 4 = 0.602. 

c. The accuracy of the current and the L/D-F schemes has a sharp rise, respectively, near 8 = 
I/V 3 =0.577 and 8 = 2/^3 = 1.155. 

The above observations can be explained by using the stability map Fig. 4.1 and Eqs. (7.58) - 
(7.60). In addition, the results shown in Fig. 7.7 also confirm the predictions given by Eqs. (7.61) 
- (7.63). 

The last test problems to be considered are those in set #8 (see Table 7.3). These problems 
share the same values of a, (i, b, K, and t, and differ in their values of At. For them, we have 8 = 
0, x = 30Ar and n = IN(r/Ar) = IN(15/x). In Fig. 7.8, the values of R M /R\ and R L /R\ are 
plotted against x. One observes that: 

a. R„IR 1 is nearly a constant (= 4) throughout the range 1 > x > 0. 

b. R l /R 1 is nearly a constant (=4) in the range 0.4 > x > 0. Its dependence on x becomes 
more and more irregular as x increases from 0.4 to 1.0. 

With the aid of Eqs. (7.14) and (7.29), observation (a) can be explained by using Eqs. (5.84) and 
(5.87). On the other hand, with the aid of Eqs. (7.14) and (7.25), observation (b) can be explained 
by using Eqs. (5.84) - (5.86) and the fact that n = IN(15/x) is relatively small when x > 0.4. 

So far in this section, no numerical results are provided for the error R 2 (J,n). This omission 
is because our discussions have focused on the comparisons of the three schemes considered, and 
there are no simple counterparts of R 2 {j,n) in the MacCormack and the L/D-F schemes. 
However, it should be noted that the errors R\(j,n), R 2 (j,n), R L (J,n), and R„(j,n) at different 
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(J,n), respectively, have been thoroughly compared against their approximations obtained by 
evaluating all the explicitly-given terms on the right sides of Eqs. (7.15), (7.22), (7.26), and 
(7.30). Without going into details, it is sufficient to state that these approximations are highly 
accurate as long as the errors they approximate are sufficiently small. 



8. FINAL REMARKS 


Using Eq. (2.2) as a model equation, the basic concepts and properties of the current 
numerical framework were described and studied in sections 1-7. By employing nontraditional 
discrete variables and taking advantage of the flexibility gained in a unified treatment of space 
and time, we were able to construct an explicit marching procedure from a single flux 
conservation principle. 

Several fundamental differences that separate the current scheme from other explicit schemes, 
and how these differences result in greater stability and accuracy for the current scheme were 
discussed and explained near the end of Section 2. Other important concepts and ideas were 
discussed in the last part of Section 1 . 

Perhaps the most intriguing results presented in the current paper are the similarities between 
the L/D-F scheme and the current scheme. It has been shown that: 

a. There is a remarkable similarity between the forms of the amplification factors of these two 
schemes. 

b. These two schemes have the same stability region on the 5-x plane. 

c. The stability condition of the current scheme, as in the case of the L/D-F scheme, is 
essentially the CFL condition and thus independent of the viscosity coefficient |i. 

d. Both schemes have no numerical diffusion in the absence of viscosity. 

e. The consistency of the current scheme, as in the case of the L/D-F scheme, requires that 
At/ Ax — » 0 as Ar, Ax — > 0. 

Since most of the above similarities are nontrivial, their existence suggests that there may be 
deeper reasons behind these similarities. Progresses made in a study along this direction will be 
reported in near future. 

Despite of the above similarities, it was shown theoretically and numerically that the current 
scheme is far superior to the L/D-F scheme in accuracy. 

In order to clarify the discussion on consistency, the concept of "strong consistency" is also 
introduced. It is shown that a scheme which is strongly consistent with the PDE being solved has 
no intrinsic advantage over a scheme which is not. 

This paper is concluded with the following odds and ends: 

a. The triangle APQR depicted in Fig. 8.1 is a boundary conservation element. Let a, (3, and y 
be constants. Let (x 0 ,to) be the coordinates of the mid-point M of QR. Let A"PQR be the 
interior of APQR. Let 
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H(x,t) = a(x-x Q ) + p (t-t 0 ) + y ( x,t ) e A"PQR 


( 8 . 1 ) 


be the approximation of u(x,t) in A"PQR. Let the flux entering APQR through the edge 
PQ and the value of u at M (denoted by u„) be given. Then a, p, and yean be determined 
by the requirements that (i) the total flux leaving APQR = 0, (ii) the flux be balanced at PQ, 
and (Hi) u(xqJ 0 ) = y = u M . Upon determining a, p, and y, one can calculate the flux 
leaving APQR through RP. 

b. In Fig. 8.2, a space-time £2 > s divided into conservation elements that are hexagons. Each 
hexagon has three incoming fluxes and three outgoing fluxes. A marching procedure can 
be defined if the outgoing fluxes are expressed in terms of the incoming fluxes. A possible 
way to do this will be described as follows. Let 

u(x,t) =? A(x-xq ) 2 + B(x-x 0 )(t-t 0 ) 

+ C(t-t 0 ) 2 + D(x-x 0 ) + £(f-f 0 ) + F (8.2) 

be the approximation of u (x,t) in the interior of a hexagoa Here A, B, C, D, E, and F are 
constants, and (x 0 ,to) are the coordinates of its center. Let 7f(.x,/) be defined by Eq. (2.4) 
and 

= 0 ( 8 . 3 ) 

Then the divergence theorem implies that the total flux leaving this hexagon = 0. Eqs. 
(8.2) and (8.3) also imply that 

B = -2a A , C = a 2 A , E = -aD + 2\aA (8.4) 

As a result, 

u(x,t) = A[(x-x 0 )-a(t-t 0 )] 2 

+ D[(x-x 0 )-a(t-t 0 )] + 2\iA(t-t 0 ) + F (8.5) 

The coefficients A, D, and F in Eq. (8.5) can be determined in terms of the three incoming 
fluxes. In turn, the outgoing fluxes can be determined by using Eqs. (2.4) and (8.5). 

c. As depicted in Fig. 8.3, a space-time can also be divided into conservation elements of 
different geometry shapes. 

d. A possible conservation element in a space-time £3 is shown in Fig. 8.4. This 
conservation element is formed by three "incoming" surfaces RSWV, QRVU, and TWVU, 
and three "outgoing" surfaces PSWT, PTUQ, and PSRQ. The "marching" direction is that 
of vf>. Because four points generally are not on a plane in £3, a surface like RSWV can be 
formed by two triangles, ARWV and ARSW. 
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In the interior of this conservation element, we may assume that 

u(x,y,t) = a(x-x 0 ) + P(j-yo) + y(f-fo) + 5 (8.6) 

where a, p, y, and 5 are constants, and (xo.yo^o) are the coordinates of the center of the 
conservation element. The parameters a, p, y, and 5 can be determined by the 
requirements that (i) the total flux leaving the conservation element = 0 and (ii) the fluxes 
be balanced at three incoming surfaces. 
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Appendix A 


In this appendix, stability will be discussed assuming 1 - x 2 + 5 * 0 and 5 > 0. 


Eq. (4.14) implies that 


Thus 


a + (0) • a_(0) 


5-(1-t 2 ) 
5 + ( 1 - x 2 ) 


(A.l) 


|a + (0)| • | o.(0) | 



An immediate result of Eq. (A.2) is 


max { |o+(0)| , | o_(0) | } > 


if x 2 > 1 and 5 > 0 
if x 2 < 1 and 5 > 0 

if x 2 > 1 and 5 > 0 


(A.2) 


(A.3) 


From Eqs. (4.33), (4.34), and (A.3), one concludes that the current scheme is unstable if x 2 > 1 
and 5 > 0. This fact coupled with a result established in subsection 4.2, i.e., the scheme is 
unstable if x 2 > 1 and 5 = 0, leads to the conclusion that the current scheme is unstable ifx 2 > 1 
and 5 > 0. 


From Eqs. (4.13) and (4.14), one has 

0 0 

o ± (0) = cos(— ) ± i | sin(— ) | if x 2 = l and 5>0 (A.4) 

Thus 

I o+(0) I = | o_(0) 1=1 if x 2 = 1 and 5 > 0 (A.5) 

Because (/) o+(0) = o_(0) is a necessary condition for the defectiveness of 2(0), (ii) in the case 
where x 2 = 1 and 5 > 0, o + (0) = o_(0) occurs only if 0 = 0, and (Hi) the matrix 2(0) is the 
identity matrix when x 2 = 1 and 5 > 0, the matrices 2(0) are nondefective for n > 0 > -n if 
x 2 = 1 and 5 > 0. As a result, Eqs. (4.33) and (A.5) imply that the current scheme is stable if 
x 2 = 1 and 5 > 0. Note that the assumption 1 - x 2 + 5 * 0 excludes the case in which x 2 = 1 and 
5 = 0. 


It was shown in subsections 4.2 and 4.3 that the current scheme is stable if either (i) x 2 < 1 
and 5 = 0, or (ii) x = 0 and 5 > 0. Thus the only stability problem left to be discussed is that in 
which 1 > x 2 > 0 and 5 > 0. 

To proceed, note that Eq. (4. 13) implies that 

h(0)] 2 + (1 -x 2 ) 2 - 5 2 = X(0) + / Y(0) 


(A.6) 



where X(0) and Y(0) are real functions defined by 


X(0) ( 1 - x 2 ) 2 [ 1 - x 2 sin 2 (-|) ] - 8 2 sin 2 (-|) 

Y(0) -2 8x(l -x 2 )sin(-|)cos(y) 


(A.7) 

(A.8) 


In the following derivations, X(0) and Y(0), respectively, may be abbreviated as X and Y. Since, 
by definition, the range of the phase angle of the principal square root is (-rc/2 , Jt/2], we have 

Vx+7Y = £ VV x 2 + Y 2 + X + i sign(K)^ X 2 + Y 2 -X j (A.9) 

where sign(Y) ^ 1 if Y > 0 and sign(K) -1 if Y < 0. By using Eqs. (4.14), (A.6) and (A.9), 
one concludes that 


a±(0) = 


5 cos(Q-) ± -jL VV x 2 + Y 2 + X 
2 V2 


+ i 


4 ignff) 


v/V x 2 + y 2 -X -t(l -x 2 )sin(-|) 


( 1 - x 2 + 5) 


It will now be shown that 

| a+(0) | < 1 and | o_(0) |<1 ifx 2 <l,8>0 and 0*0 


(A. 10) 


(A. 11) 


Proof. Using Eq. (A. 10), one has 

[| O + ( 0 )| 2 + | O _( 0)| 2 ] ( 1 - t 2 + 5) 2 

= 2 [ 8 2 cos 2 (y ) + x 2 ( 1 - x 2 ) 2 sin 2 ( y) + ^X 2 + Y 2 ] (A. 12) 

Combining Eqs. (A.l) and (A. 12), one obtains 

[1- |a + (0)| 2 ] [1- |a_(0)| 2 ](l-x 2 + 8) 2 

= 2 { ( 1 - x 2 ) 2 [ 1 - x 2 sin 2 (y) ] + 8 2 sin 2 (|-) - Vx 2 + Y 2 } (A. 13) 

Because (i) 

{(1— x 2 ) 2 [1— X 2 sin 2 (y ) ] + 8 2 sin 2 (y ) - ^X 2 + Y 2 } 
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(A. 14) 


X { ( 1 -x 2 ) 2 [ 1 -t 2 sin 2 (y ) ] + 5 2 sin 2 (-|) + Vx 2 + Y 2 } 
= 4 6 2 ( l - t 2 ) 3 sin 2 (y) 


(l-x 2 ) 2 [l-t 2 s in 2 (|)] + 5 2 sin 2 (|) + T/x T TY r > 0 if x 2 < 1 (A.I5) 

and (Hi) 

45 2 (1 -x 2 ) 3 sin 2 (|-) > 0 if x 2 < 1 , 5>0 and 9*0 , (A.16) 

one concludes that the expression on the the right side of Eq. (A. 13) > 0 if x 2 < 1, 5 > 0 and 
0 5*0. Thus 


[1 - Ja + (0)j 2 ] [ 1 - |a_(8)| 2 ] >0 if x 2 < 1 , 5 > 0 and 0*0 (A.17) 

According to Eq. (A.2), |a + (0)j • |a_(0)| <1 if x 2 < 1 and 5>0. This inequality combined 
withEq. (A.17) implies Eq. (A. 11). Q.E.D. 

Note that (i) 

I o+ (0)| = 1 and |o_(0)| = l ~g~~7T~Tv l < 1 if x 2 < 1 and 5 > 0 (A.18) 

and (ii)Q(O) is diagonal and thus nondefective. According to Eqs. (4.33) and (4.34), (i), (ii), and 
Eq. (A. 11) imply that the current scheme is stable if t 2 < 1 and 5 > 0. Combining this and other 
results obtained earlier, one concludes that, assuming 8^0 and 1 — t 2 + 8^0, the current scheme 
is stable if and only if x 2 < 1 . 
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Appendix B 


The main L/D-F scheme is defined in Eq. (3.19). It will be supplemented by the starting 
scheme 


±A 

At 


-2.g . 0 

" 1 (Ax) 2 


2 Ax 


(B.l) 


Because the moving mesh shown in Fig. 2.1(a) is used in the current discussion, coefficient a is 
replaced by a - b in Eq. (B.l). We will also assume periodic conditions: 


u ] ~ U j+K (y -0, ±1, ±2, • • • , n =0, 1, 2, • • ■ ) 
As a result of Eq. (B.2), 

u] = £ ul 

* = 0 


where are defined by Eq. (4.2) and 

u n k ^ s' «7 W 

I =0 


(B.2) 


(B.3) 


(B.4) 


Substituting Eq. (B.3) into Eq. (3.19) and using Eqs. (2.17), (2.18), and (4.7), one obtains 
O + y ) «lfc +1 + (2i T Sin0* - 5cos0*) - ( 1 - y ) m* _I = 0 , 11 = 1,2.3, ••• (B.5) 

Note that Eq. (B.5) can be expressed as 

iy +1 = M(Q k )tk , n=0, 1,2, ••• (B.6) 


where 


and 




r «z +1 


u™ 

Hk 


M(0) = 


5cos6-2/xsin6 
1 + (5/2) 

1 


1 -(5/2) , 
1 +( 8 / 2 ) 

0 


(B-7) 


(B.8) 


As a result of Eq. (B.6), M(Q k ), k = 0, 1, 2, • • • , AT— 1, can be referred to as the amplification 
matrices of the L/D-F scheme. Also we have 

S’ = [M(Q k )] n tk . «=0, 1,2, ••• (B.9) 
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The eigenvalues of M(0), i.e., 


(y )cos0 - i x sin0 ± '\J (y-cos0 - ix sin0) 2 + 1 - ( — ) 2 

; (B.10) 

1 + f 

will be referred to as the amplification factors of the L/D-F scheme. 

Because Af(0) is not diagonal, it is defective if and only if A + (0) = /4_(0). With the aid of this 
fact and Jordan’s theorem [p.362, 6], Eqs. (B.7) and (B.9) can be used to show that 

f hk+[A + (Qk)f +h k -[A_(Q k )] n if i4 + (0*)^i4_(0 t ) 

^ = 1 , (B.l 1) 

l h k+ [ A + (0*)]- + nh k .[ /U0*))"- 1 if A + (9 k ) = A.(9 k ) 

where h k + and h k _ are constants to be determined by ~d k . In this paper, the L/D-F scheme is said 
to be stable if and only if, for any K > 3 and any specification of h k+ and /i*_, k = 0, 1, 2, • • • , 
tf-1 (i.e., any specification of u'j and «],;' = 0, 1, 2, • • • , K- 1), u ], j = 0, ± 1 , ± 2 , • • • , remain 
bounded as n -» +« with the parameters x and 5 being held constant. From Eq. (B.l 1), one can 
conclude that the L/D-F scheme is stable if and only if, for any K >3 and any k = 0, 1, 2, 

AT — 1 , we have 

max { |i4 + (0*)| , |/4_(0*)| } < 1 if >l + (0*) ^y4_(0*) (B.12) 

and 

M+(9*)| < 1 if A + (9*) = >1 _(©*) (B.l 3) 

In the following derivations, we shall prove this: In the case where 8 > 0, the L/D-F scheme 
is stable if and only if x 2 < 1 . In the case where 5 = 0, it is stable if and only if x 2 < 1 . 


Proof : It is easy to show that 
max { |i4 + (-y)| , |A_(-|)| } 
and 


1x1 +Vx 2 -1 +(5/2) 2 

1 + ( 8 / 2 ) > 


if x 2 > 1 and 8>0 (B.14) 


^+(j) =^-(y) ^ M+(y)| = = 1 if x 2 = l and 8 = 0 (B.15) 

Because 0* = n/2 if K = 4, 8, 12, • • • , and k =K/ 4, Eqs. (B.12) - (B.15) imply that the L/D-F 
scheme is unstable if either (i) x 2 > 1 and 8 > 0, or (ii) x 2 > 1 and 8 = 0. 


Next we observe that 



|A + (9)| = | A _(9)| = 1 


and 


[ if x 2 < 1 and 5 = 0 (B.16) 


A + (Q) * A_(0) 

It follows from Eq. (B. 12) that the L/D-F scheme is stable if x 2 < 1 and 5 = 0. 

The proof is completed by showing that the L/D-F scheme is stable if x 2 < 1 and 5 > 0. To 
proceed, from Eq. (B.10), we obtain 


and 


A + (0)^_(0) 


5/2-1 
5/2 + 1 


A±(Q) = 


+ 


4 COS0 ± -pr VV X 2 +Y 2 +X 

2 ^2 j 

±-L sign(F) VV x 2 + Y 2 -T - xsin0 
V2 


(1 + 6 / 2 ) 


(B.17) 


(B.18) 


where X and Y, respectively, are the abbreviations of 

X(0) (1 -x 2 sin 2 0) - (|-) 2 sin 2 0 , and Y(0) ^ -5xsin0cos0 (B.19) 

Also sign(T) 1 if Y > 0 and sign(F) d - -1 if Y <0. Eq. (B. 18) implies that 

[|A + (0)| 2 + |A_(0)| 2 ] (1 + -|) 2 = 2[(|-) 2 cos 2 0 + x 2 sin 2 0 + Vx 2 + Y 2 ] (B.20) 


By using Eqs. (B.17) and (B.20), one concludes that 

[1- |A + (0)| 2 ][1- |a_(0)| 2 ](i + |) 2 

= 2 [ 1 - x 2 sin 2 0 + ( ) 2 sin 2 0 - Vx 2 + Y 2 ] (B.21) 

Note that 

[ 1 -x 2 sin 2 0 + ( f sin 2 0 - ^X 2 + Y 2 ] [ 1 - x 2 sin 2 0 + ( f sin 2 0 + Vx 2 + Y 2 ] 

= 6 2 (1 -x 2 )sin 2 0 (B.22) 


and 
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if x 2 £ 1 and 8 > 0 (B.23) 


[ 1 - 1 2 sin 2 0 + (|) 2 sin 2 0 + Vx 2 + y 2 ] > 0 

Combining Eqs. (B.21) - (B.23), one concludes that 

[ 1 ~ M+(0)| 2 ] [ 1 - |A_(9)| 2 ] = 0 if x 2 = 1 and 5 > 0 
and 


(B.24) 


[ 1 - |A+(0)| 2 ] [ 1 - |A_(0)| 2 ] > 0 if x 2 < 1 , 5>0 and sin0*O (B.25) 
An immediate result of Eq. (B.17) is 

|A + (0)| • |A_(0)| < 1 if 5 > 0 
From Eqs. (B.24) - (B.26), one concludes that (i) 
max { |j4 + (0)| , |A_(0)| } = 1 


(B.26) 


and 


)• if x 2 = 1 and 8 > 0 (B.27) 


M+(0)| * |A_(0)| 

and (ii) 

M+(0)| < 1 and |A_(0)| < 1 if x 2 < 1 , 8 > 0 and sin0*O (B.28) 
Furthermore, Eq. (B.10) implies that 

max { |A + (0)| , |A_(0)| ) = 1 


and 


|A + (0)| * |A_(0)| 


f if sin0 = O and 5 > 0 (B.29) 


By using Eqs. (B.12), (B.13), and (B.27) - (B.29), one reaches the conclusion that the L/D-F 
scheme is stable if x 2 < 1 and 8 > 0. Q.E.D. 

Next we shall study the accuracy of the numerical solutions obtained by using the L/D-F 
scheme. Let 


., ^+(0) - 1 + — (l-cos0) + txsinO 

h(Q ) ? 

A + (Q) - A_(0) 


(A + (0)*A_(0)) (B.30) 


and 


lb.(j,n,k) =f b k e‘ jQk { [ 1 - A(0 t )] [A + (0 t )] B + A(0*)[A_(0*)]" } 
(A + (0 t ) * A_(0*)) 


(B.31) 
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where b k are defined in Eq. (5.8). With the aid of Eqs. (B.3) and (B.l 1), and the assumption that 


^ + (0*) ^ A_(0 a ) , k=0, 1, 2, .... AT— 1 


(B.32) 


it can be shown that the solution to the finite-difference problem defined by Eqs. (3.19), (B.l), 
and (B.2) is given by 

def * _1 

u] = &.(/>«) - £ «*.(/>.*) (B.33) 

*= o 


Note that 


A + (0) = 1 


and 


A_(0) = - 


1-5/2 
1 + 5/2 


(B.34) 


Thus there is a neighborhood of 0 = 0 on the complex 0 -plane in which A + (0) * A_(0) and thus 
A(0) is defined. 


Combining Eqs. (5.53) and (B.31), one has 

Ml (/.«>&) - u a (J,n,k) = b k e lj9k A L (n,Q k ) (B.35) 

where 

a,(«,0) ^ [i-/t(0)nu + (0)r-M a (0)n + m[[A.(Q)] n -[A a (Q)r) (b.36) 

In the following, we shall study A L (n, 0) in the neighborhood of 0 = 0. 

Eqs. (B.10) and (B.30) imply that 

h(d) = l(i-^-)|_ x 2e 2 + -^-(i-2t 2 )0 3 

+ [ 4 x 2 ( 4 - 9 x 2 ) + 3 S 2 ( 1 5 x 4 - 12 x 2 + 1 ) ] 0 4 + 0(0 5 )j (B.37) 

Eq. (B.37) is reduced to 

2 2 2 

h(Q) = ~^-0 2 + ~ (4 4 ^ 9T ) Q 4 + 0(Q 5 ) if 5 = 0 (B.38) 

and 

A(0) = “5 2 0 --^-)0 4 + 0(0 5 ) if x = 0 (B.39) 

According to Eq. (5.12), the sum of the upper elements in 7/ + (0) and 7/_(0) is 1. Thus the role 
played by the upper element of 7?_(0) in Eq. (5.21) is similar to that played by /t(0) in Eq. (B.36). 
An inspection of Eqs. (5.32) - (5.34) and (B.37) - (B.39) reveals that the upper element in /?_(0) 
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is (i) smaller than /t(0) by one order of 0 if x 5 * 0, or (ii) smaller than h(Q) by two orders of 0 if 
5 = 0, or ( Hi) in the same order of 0 with h (0) if x = 0. Thus, if x ^ 0, the influence of the spurious 
part is less noticeable in the current scheme than in the L/D-F scheme. Note that the form of h(Q) 
given in Eq. (B.30) is dependent on the forms of both the main scheme Eq. (3. 19) and the starting 
scheme Eq. (B.l). It is an accident of this combined influence of the main and starting schemes 
that /i(0) and the upper element in 7/_(0) are in the same order of 0 if x = 0, even though du /dt is 
approximated by a one-sided difference formula in Eq. (B.l). 

Next, with the aid of Eqs. (5.17) and (B.10), it can be shown that 


«f<9> * 4 1 ? - 1 = 


A,® 


3x 


+ 96 l d-t 2 )[^-( 5 5 2 - 25 + 4) + 2(l-f S 2 )] + 35x 2 j-0 4 + 0(0 5 )(B.4O) 


V 


J 


and 


e,-(0) ^ 


A -(d) 


-( 


1-5/2 

1+5/2 


K<6) 


- 1 = 2/x0 + [-^(2-x 2 ) 
4 


2x 2 ] 0 2 + 0(0 3 ) (B.41) 


Obviously, the roles played by e i+ (0) and £^-(0) in the L/D-F scheme, respectively, are similar to 
those played by e+(0) and e_(0) in the current scheme. A comparison between Eqs. (5.36) and 
(B.40) reveals that, [a+(0)] 2 approximates A a (Q) to the second order in 0 while A+(Q) 
approximates A a (0) only to the first order in 0. Note that the amplification factors are completely 
determined by the main marching scheme. Thus, the above difference in accuracy has nothing to 
do with how the starting scheme is defined. 


By using the arguments that were used to establish Eqs. (5.40) and (5.45), it can be shown 
that (i) 


[A + (0)r-[A fl (0)] n + n [ A a (0) ]" e t+ (0) (B.42) 

if n | £/.+(©) | <c 1, and (ii) , 

[A_(e)} n -[A a (0)f = -[A a (0)]" (B.43) 

if 5 > 0, 1 0 1 is sufficiently small, and n is sufficiently large. 

Let 


r L +(n,k) 


dg [l-h(B)]{[A + m" - [A a (0)p 

[A a m n 


(B.44) 
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n I 



and 


r L -(n,k) 


def 


m[[Ajfi)r - [4,(0)n 

[A a m n 


(B.45) 


Then Eqs. (5.61) - (5.63) follow directly from Eqs. (5.53), (5.60), (B.35) - (B.37), (B.40), (B.42), 
and (B.43). 


For the special case in which 5 = 0, Eq. (B.43) is not applicable. However, since |.A_(0)| = 
I Ai(9) | = 1 if x 2 < 1 and 5 = 0, one has 


M-(9)]" ~ [i4 a (9)r 

[A a m n 


< 2 


if x 2 < 1 and 5 = 0 


(B.46) 


With the aid of Eqs. (B.38) and (B.45), Eq. (B.46), in turn, can be used to obtain Eq. (5.85). 
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Appendix C 

Eqs. (5.26) - (5.28) and the assumption «|e + (9)| <e 1 of Eq. (5.40) will be discussed in this 
appendix. 

Assuming Eqs. (5.23) and (5.24), first we shall show that 

a. ©(9) * 0 if 7 C > 1 0 1 ; and 

b. co( 7 i) = 0 if and only if 1 - x 2 - 8 = 0 and 1 > x > 0 


Proof. By using Eqs. (4. 1 3), (4.14), (5.23), and (5.24), it is easy to show that 


<m 


2(1 -X 2 ) 

-x 2 + 5 


(C.l) 


Let I be the 2 x 2 identity matrix. Then the determinant of the matrix [ 2(9) - a_(0) I ] vanishes. 
This fact coupled with Eqs. (4. 12), (5.23), (5.24), and (C.l) implies that 


to(0) ^ 0 if 7 c > 0 > — k and 1 — x 2 — 8 ^ 0 (C.2) 

Let 1 -x 2 - 8 = 0. Then 5 = 1 - x 2 > 0 . This inequality combined with Eq. (4.13) implies that 
the real part of rj(0) is positive if it > [ 0 1 . As a result, the principal square root ^[r|(0 )] 2 = rj( 0 ) 
if 7 i > j 0 1 . Combining Eqs. (4.14) and (5.24), one can conclude that o_( 0 ) = 0, and 

©(0) = cos(0/2) — i xsin(0/2) * 0 , if it > | 0 | and 1— x 2 -S = 0 (C. 3 ) 

Statement (a) ia a result of Eqs. (C.2) and (C.3). Statement (b) follows directly from Eq. 
(C.2) and the fact that 


©00 - («/2)( | x [ -x) if 1— x 2 — 8 = 0 (C.4) 

Q.E.D. 

Next we shall prove Eqs. (5.26) - (5.28) for any 0 with jc > 0 > -it and o + (0) * o_(0). To 
proceed, note that a + (0) + a_(9) = the trace of 2(9). By using Eqs. (4.12) and (5.24), one obtains 

1 5 

= ( 73 ^ 27 ^ )[ c °s(9 /2 ) + txsin(0 / 2)] + a + (e) ( re > 0 > -re ) (C.5) 

Because c + ( 0 ) * a_( 0 ), ^+( 0 ) and ^1(0), respectively, are the eigenvectors of 2 ( 0 ) with 
eigenvalues o + (0) and o_(0). With the aid of Eqs. (4.12), (4.22), (5.24), (5.25), and (C.5), it can 
be shown that 


and 


§ 21 (9) = «^i(9)g 11 (0) 


(C. 6 ) 
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T 1 I 


gl2(9) “ *§2(6) #22(9) 


(C.7) 


By choosing g n (0) = £ 22 ( 8 ) = 1- Eqs. (4.21), (C. 6 ), and (C.7) imply that 


G(0) = 


1 / § 2 ( 9 ) ' 

i§ i(0) 1 


(C. 8 ) 


Because ^+(0) and^_(0) must be linearly independent (i.e., G(0) must be nonsingular) if o+(0) * 
a_(0), Eqs. (5.26) - (5.28) follow from Eqs. (C. 8 ) and (5.1 1) immediately. 

As a preliminary for a discussion on the assumption n |e+(0)| <*: 1 of Eq. (5.40), we shall 
establish the following inequality: 


Lemma. Let x be a complex number. Let n > 0 be an integer. Then 

|(1 + x) 1 - 1 1 <*: 1 , / = 1, 2, 3, .... n 


if and only if 


n \x\ 1 


(C.9) 


(C.10) 


Proof. Let the real numbers r, <J>, p, and \| / be defined by the conditions: 

x = re 1 *, r> 0, rc > <|> > -7C (C.ll) 

and 

l+* = p<; iv , p > 0 , jt>y>-JC (C.12) 

(See Fig. C.l). It follows from Eq. (C.12) that 

| ( 1 + x - 1 1 = ^ ( p 7 - 1 ) 2 + 2 p 1 ( 1 - cos(fy) ) 

By using Eq. (C.l 3) and the fact that 1 - cos (lx\f) > 0, one concludes that Eq 
only if 

(p ; - 1 ) 2 <c 1 , / = 1, 2, 3, .... n 

and 

1 - cos(/\|/) <£ 1 , / = 1, 2, 3, .... n (C.15) 

Note that Eq. (C.14) implies that p' is very close to 1 for / = 1, 2, 3, .... n. As a result, the second 
term under the radical sign on the right side of Eq. (C.l 3) will be small compared with 1 for 
l = 1, 2, 3 n if and only if Eq. (C.15) is true. Because rc > y > -jc, Eq. (C.15) is equivalent to 


(C.l 3) 

. (C.9) is true if and 
(C.14) 
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Also, Eq. (C. 14) is equivalent to 


A result of Eq. (C.17) is 


« |>|/| <ec 1 


| e | « 1 where e ^ p" _ \ 


(C.16) 


(C.17) 


P = ( 1 + e) " = 1 + -^ 

n 


(C.18) 


With the above preparations, first we shall show that Eq. (C.10) is a result of Eq. (C9) 
According to Fig. C.l, 


(C.19) 


1*1 - r = ^ p 2 + 1 -2pcosxjr 
By using Eqs. (C.16) - (C.19), one concludes that 

»W ±»V + l-2p(l-^/2) =„V (p _ 1) 2 +pv 2 ± n'T(e/n) 1 +(l + e/n) v 2 
= ^ e 2 + (1 +e/«)(«vjr) 2 <*; i 
i.e., Eq. (C.10) is true. 

Next, assuming Eq. (C.10), we shall prove Eq. (C.9) by induction. Obviously |(1 +x) 1 - 1 
= \x | c 1 if /= 1. Let 7 0 be an integer such that n > 7 0 > 1 and 


|(l+*y-l| <e 1 , 7=1,2, 3, .... / 0 

As a result of Eq. (C.20), we have 

K 1+ *)'l < 2 , l=l, 2,3 / 0 

With the aid of Eqs. (C.10) and (C.21), we have 

K i +*) /o+1 - 1 1 s i* s o+*)'i ^ 1*1 z \i+x\ l < (2/ 0 

1=0 1=0 

Q.E.D. 

With the aid of the above lemma, we now show that 

f <7+(0) J 27 = [A a (0)] / , 7 = 1,2 n („> 0) 

if and only if 

n le^O)! « 1 


(C.20) 


(C.21) 


(C.22) 


(C.23) 
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Proof: By replacing n with / in Eq. (5.38) and recalling the definition of the sign "=" given in 
Section 5, it is seen that Eq. (C.22) is true if and only if 

| [ 1 + e + (0) ] / - 1 1 « 1 , 1 = 1, 2, 3 n (C.24) 

According to the above lemma, Eq. (C.24) is true if and only if Eq. (C.23) is true. Q.E.D. 
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Appendix D 


A version of the MacCormack scheme forEq. (2.2) [p.163, 3] is 


Predictor: 

Corrector: 


+1 " ~ fa ** ( w ? +1 ~ ) + ( «"+i ~2u]+ «;_! ) 


(D.l) 


«r* = + u] + ' - 


(a-b)At f 


( u] +l - Uj+} ) 


Ax 


* +«?-/>] 


(D.2) 


Because the moving mesh depicted in Fig. 2.1(a) is used in the current discussion, coefficient a is 
replaced by a-b in Eqs. (D. 1) and (D.2). In the above version, a forward difference is employed 
m the predictor step for du/dx and a backward difference is used in the corrector step. The 
alternate version employs a backward difference in the predictor step and a forward difference in 
the corrector step. For both versions, the predictor and corrector steps can be combined to yield 

Uj = J (j +x ) u J-2 + +x + t2 -~Y~-^-) u ]-i + ( 1 -|’ + ^-5 2 -x 2 ) m ; 


4(f^ 2 + T-fK^!(7-^ 


(D.3) 


Let 


and 


>1(9) « 0-j + 1 |-s’-T 2 ) + (f + T 3 --^)co S e 

“ i'T(l--)sin0 + -y— cos20 - -^sin20 
-6 lo 4 


u M (J,n,k) b k e i} ^[A{Q k )T 


(D.4) 


(D.5) 


where b k are defined in Eq. (5.8). Then the solution to the finite difference problem defined by 
Eqs. (D.3) and (B.2) is 


u j = &*(/.") ^ Z «»(/.«.*) 

k = 0 


*■-1 


(D.6) 


Combining Eqs. (5.17) and (D.4), it can be shown that 
£(0) ~ A a (0) “ 1 = - T(1 6 ~ - 0 3 + ^-[8(1-6x 2 )-6t 2 (1-x 2 )]0 4 + 0(6 5 ) (D.7) 


By using the arguments used to establish Eq. (5.40), one can show that 
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M(9)l" - [A a (Q)] n = «Ma(e)] B e(0) if n |e(0)| c 1 


(D.8) 


Eq. (5.65) follows directly from Eqs. (5.64), (D.5), (5.53), (D.8), and (D.7). 

We conclude this appendix with a discussion on the operation count of the MacCormack 
scheme. Because the coefficients in front of the mesh variables on the right side of Eq. (D.3) are 
constants, they need not be reevaluated during the numerical marching. Thus, for each j, the 
MacCormack scheme requires 5 multiplications and 4 additions to advance one time step. 
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Appendix E 


The proofs for several assertions made in Section 6 are provided here. 

We begin with the assertion: Let p > 0. Then the Lax stability of the MacCormack scheme 

for solving Eq. (2.2) requires that the mesh be refined such that the parameter 5 remains bounded 
as At, Ax — » 0. 


Proof. Because (i) the sum of the terms involving 8 2 on the right side of Eq. (D.4) = 
(cos0 - 1 ) 2 5 2 /8, and (ii) ^im Q x/5 = 0, one concludes that, for any 0 with cos0*l, the 

amplification factor A (0) of the MacCormack scheme will become unbounded as At, Ax 0 if 5 
becomes unbounded as At, Ax — » 0. Since uniform-boundedness of the spectral radius of the 
amplification matrix is necessary for the Lax stability [p.70, 4], the proof is completed. Q.E.D. 

Next we show that Eqs. (6.29) and (6.30) are true if u = u(x,t) and v = v(x,t) satisfy Eq. 
(6.27), i.e., w^x.r) = 0. To proceed, note that 


1 > 


1 -x 2 - 5 
1 - x 2 + 8 

follow directly from Eq. (6.18). Also, 


| and 1 > 1 . T — > 0 

1 1 - x 2 + 5 


x± 


8(1 tx) 

1 - x 2 + 8 

8( 1 Tx) 
1 -x 2 + 8 


1 i X 


At 


xO[( Ax) 4 ] = 


a ( Ax ± a At) ± 4 p 


1 — x 2 


1 - x 2 + 8 


(E.1) 


x O [ (Ax) 2 ] (E.2) 


x± 


1 - x 2 - 8 

At 


x O [ (Ax) 4 ] 


= \ a [ (A*) 2 - a 2 (Ar) 2 - 4 p (Ar) ] ± 4p(Ax Ta At) — tz 5 

L 1 - x 2 + 8 


xtf[(Ax)] 


(E.3) 


With the aid of Eqs. (6.23) and (E.l) - (E.3), and the fact that dw/dx = 0 if w(x,t) = 0, one 
concludes that Eq. (6.29) is true if u = u(x,t) and v = v(x,t) satisfy Eq. (6.27). 

Next, we have 


0 and Ax 8 


0 


as At /Ax 


0 


Also, 


(1 -x 2 +8) 2 (Ax) 2 = (Ax-x 2 Ax + 8Ax) 2 


— — T g + 5)2 x O [ (Ar) 2 ] = (A * + 8 Ax ) 


4p 


xO[(Ar)] 


(E.4) 

(E.5) 

(E.6) 
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"T Ii 


i±ja-i 2 > 


= , 1 ^( 1 -#) 

4M- 


( 1 - t 2 - 5 ) x <9 [ (Ax) 2 ] = O[(At) 2 ] + <9[(Af) 2 ]-5A*x0[(A*)] 


(E.7) 

(E.8) 


Eq. (6.30) now follows directly from Eqs. (6.24) and (E.4) - (E.8). 

If and the rule of mesh refinement is such that 5 remains bounded as At, Ax — » 0, then 
(i) At = O [(Ax) 2 ] and (ii) At /Ax 0 as At, Ax -» 0. With the aid of Eqs. (E.l) - (E.8), it is easy 
to show that Eq. (6.31) is true if u = u(x,t) and v = v(x,t) satisfy Eq. (6.27) and the rule of mesh 
refinement is such that 8 remains bounded as Ar, At — > 0. 
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a 

D 

b 

K 

t 

Re 

x 0 

x* 

% 

# 1 

1.0 

0.1 

0 . 

30 

0.5 

1/12 

0.048002 

0.048011 

0.1064 

#2 

1.0 

0.1 

0 . 

30 

1.0 

1/12 

0.048002 

0.018026 

0.1064 

#3 

1.0 

0.1 

0 . 

60 

0.5 

1/24 

0.024042 

0.024043 

0.05364 

#4 

1.0 

0.01 

0 . 

30 

3.1 

5/6 

0.40299 

0.40309 

0.6380 


Table 7.1. Definitions of the problem sets #1 - #4 and the corresponding 
values of Re, x 0 , x*, and To- 
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11 

K 

n 

At 

5 

x+ 

X* 


#5 

1.0 

0.01 

30 

390 

l/(45^3) 

O.8/V3 

0.4472 

0.4468 

0.8014 

#6 

1.0 

0.01 

30 

312 

l/(36^3) 

I/V3 

0. 

0. 

0.7435 


Table 7.2. Definitions of the problem sets #5 and #6, and the corresponding 

values of 8, x + , x*. and x + . 
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0 . 
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15. 

#8 

1.0 

0. 

0 . 

30 

0.5 


Table 7.3. Definitions of problem sets #7 and #8. 
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C our ant number, r 

Figure 7.2— Accuracy of test problems in Set #2 (t - 30 At). 
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Rgure 7.3 — Accuracy of test problems In Set #3 (t - 60 At). 
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Figure 7.5— Accuracy of test problems In Set #5 (t - 2(1-b)/(3V3)). 
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Figure 7.8 — Accuracy of test problems in Set #8 (t > 30At). 
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Figure 8.1 —Boundary conservation element (APQR). 


Figure 8.2— Hexagons as conservation elements. 



Figure 8.3.— Space-time E 2 divided into conservation elements of 
different geometric shapes. 



Figure 8.4 —Conservation element in space-time E 


Imaginary axis 



Real axis 


Figure C.1 —Geometric relations among the parameters r, 4>, p 
and i|j defined by Eqs. (C.1 1 ) and (C. 12). 
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